diff --git a/docs/Manifest.toml b/docs/Manifest.toml index ae86df997..210f1e18e 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -1,9 +1,24 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.11.1" +julia_version = "1.11.5" manifest_format = "2.0" project_hash = "e4d965798d55f5903483b71f545533e62ea32fdf" +[[deps.ADTypes]] +git-tree-sha1 = "be7ae030256b8ef14a441726c4c37766b90b93a3" +uuid = "47edcb42-4c32-4615-8424-f2b9edc5f35b" +version = "1.15.0" + + [deps.ADTypes.extensions] + ADTypesChainRulesCoreExt = "ChainRulesCore" + ADTypesConstructionBaseExt = "ConstructionBase" + ADTypesEnzymeCoreExt = "EnzymeCore" + + [deps.ADTypes.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" + EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869" + [[deps.AMD]] deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse_jll"] git-tree-sha1 = "45a1272e3f809d36431e57ab22703c6896b8908f" @@ -28,9 +43,9 @@ weakdeps = ["ChainRulesCore", "Test"] [[deps.AbstractGPs]] deps = ["ChainRulesCore", "Distributions", "FillArrays", "IrrationalConstants", "KernelFunctions", "LinearAlgebra", "PDMats", "Random", "RecipesBase", "Reexport", "Statistics", "StatsBase", "Test"] -git-tree-sha1 = "010e6b5e06a9a6b1be6b146a7f29d5848490cc64" +git-tree-sha1 = "8a05cefb7c891378c89576bd4865f34d010c9ece" uuid = "99985d1d-32ba-4be9-9821-2ec096f28918" -version = "0.5.21" +version = "0.5.24" [[deps.AbstractMCMC]] deps = ["BangBang", "ConsoleProgressMonitor", "Distributed", "LogDensityProblems", "Logging", "LoggingExtras", "ProgressLogging", "Random", "StatsBase", "TerminalLoggers", "Transducers"] @@ -45,12 +60,13 @@ version = "0.4.5" [[deps.Adapt]] deps = ["LinearAlgebra", "Requires"] -git-tree-sha1 = "50c3c56a52972d78e8be9fd135bfb91c9574c140" +git-tree-sha1 = "f7817e2e585aa6d924fd714df1e2a84be7896c60" uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -version = "4.1.1" -weakdeps = ["StaticArrays"] +version = "4.3.0" +weakdeps = ["SparseArrays", "StaticArrays"] [deps.Adapt.extensions] + AdaptSparseArraysExt = "SparseArrays" AdaptStaticArraysExt = "StaticArrays" [[deps.AdvancedMH]] @@ -77,9 +93,9 @@ uuid = "66dad0bd-aa9a-41b7-9441-69ab47430ed8" version = "1.1.3" [[deps.ArgCheck]] -git-tree-sha1 = "680b3b8759bd4c54052ada14e52355ab69e07876" +git-tree-sha1 = "f9e9a66c9b7be1ad7372bbd9b062d9230c30c5ce" uuid = "dce04be8-c92d-5529-be00-80e4d2c0e197" -version = "2.4.0" +version = "2.5.0" [[deps.ArgTools]] uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" @@ -99,9 +115,9 @@ version = "3.5.1+1" [[deps.ArrayInterface]] deps = ["Adapt", "LinearAlgebra"] -git-tree-sha1 = "017fcb757f8e921fb44ee063a7aafe5f89b86dd1" +git-tree-sha1 = "9606d7832795cbef89e06a550475be300364a8aa" uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -version = "7.18.0" +version = "7.19.0" [deps.ArrayInterface.extensions] ArrayInterfaceBandedMatricesExt = "BandedMatrices" @@ -187,10 +203,10 @@ uuid = "62783981-4cbd-42fc-bca8-16325de8dc4b" version = "0.1.6" [[deps.Bzip2_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "8873e196c2eb87962a2048b3b8e08946535864a1" +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "1b96ea4a01afe0ea4090c5c8039690672dd13f2e" uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0" -version = "1.0.8+4" +version = "1.0.9+0" [[deps.CPUSummary]] deps = ["CpuId", "IfElse", "PrecompileTools", "Static"] @@ -200,21 +216,21 @@ version = "0.2.6" [[deps.Cairo_jll]] deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"] -git-tree-sha1 = "009060c9a6168704143100f36ab08f06c2af4642" +git-tree-sha1 = "fde3bf89aead2e723284a8ff9cdf5b551ed700e8" uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a" -version = "1.18.2+1" +version = "1.18.5+0" [[deps.CalibrateEmulateSample]] deps = ["AbstractGPs", "AbstractMCMC", "AdvancedMH", "ChunkSplitters", "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.6.1" +version = "0.7.0" [[deps.ChainRulesCore]] deps = ["Compat", "LinearAlgebra"] -git-tree-sha1 = "1713c74e00545bfe14605d2a2be1712de8fbcb58" +git-tree-sha1 = "06ee8d1aa558d2833aa799f6f0b31b30cada405f" uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "1.25.1" +version = "1.25.2" weakdeps = ["SparseArrays"] [deps.ChainRulesCore.extensions] @@ -239,9 +255,9 @@ version = "0.8.5" [[deps.CodecZlib]] deps = ["TranscodingStreams", "Zlib_jll"] -git-tree-sha1 = "bce6804e5e6044c6daab27bb533d1295e4a2e759" +git-tree-sha1 = "962834c22b66e32aa10f7611c08c8ca4e20749a9" uuid = "944b1d66-785c-5afd-91f1-9de20f533193" -version = "0.7.6" +version = "0.7.8" [[deps.CommonSubexpressions]] deps = ["MacroTools"] @@ -256,9 +272,9 @@ version = "1.0.0" [[deps.Compat]] deps = ["TOML", "UUIDs"] -git-tree-sha1 = "8ae8d32e09f0dcf42a36b90d4e17f5dd2e4c4215" +git-tree-sha1 = "3a3dfb30697e96a440e4149c8c51bf32f818c0f3" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "4.16.0" +version = "4.17.0" weakdeps = ["Dates", "LinearAlgebra"] [deps.Compat.extensions] @@ -291,9 +307,9 @@ uuid = "88cd18e8-d9cc-4ea6-8889-5259c0d15c8b" version = "0.1.2" [[deps.ConstructionBase]] -git-tree-sha1 = "76219f1ed5771adbb096743bff43fb5fdd4c1157" +git-tree-sha1 = "b4b092499347b18a015186eae3042f72267106cb" uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" -version = "1.5.8" +version = "1.6.0" weakdeps = ["IntervalSets", "LinearAlgebra", "StaticArrays"] [deps.ConstructionBase.extensions] @@ -331,9 +347,9 @@ version = "1.7.0" [[deps.DataStructures]] deps = ["Compat", "InteractiveUtils", "OrderedCollections"] -git-tree-sha1 = "1d0a14036acb104d9e89698bd408f63ab58cdc82" +git-tree-sha1 = "4e1fe97fdaed23e9dc21d4d664bea76b65fc50a0" uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -version = "0.18.20" +version = "0.18.22" [[deps.DataValueInterfaces]] git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6" @@ -362,6 +378,56 @@ git-tree-sha1 = "23163d55f885173722d1e4cf0f6110cdbaf7e272" uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" version = "1.15.1" +[[deps.DifferentiationInterface]] +deps = ["ADTypes", "LinearAlgebra"] +git-tree-sha1 = "210933c93f39f832d92f9efbbe69a49c453db36d" +uuid = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63" +version = "0.7.1" + + [deps.DifferentiationInterface.extensions] + DifferentiationInterfaceChainRulesCoreExt = "ChainRulesCore" + DifferentiationInterfaceDiffractorExt = "Diffractor" + DifferentiationInterfaceEnzymeExt = ["EnzymeCore", "Enzyme"] + DifferentiationInterfaceFastDifferentiationExt = "FastDifferentiation" + DifferentiationInterfaceFiniteDiffExt = "FiniteDiff" + DifferentiationInterfaceFiniteDifferencesExt = "FiniteDifferences" + DifferentiationInterfaceForwardDiffExt = ["ForwardDiff", "DiffResults"] + DifferentiationInterfaceGPUArraysCoreExt = "GPUArraysCore" + DifferentiationInterfaceGTPSAExt = "GTPSA" + DifferentiationInterfaceMooncakeExt = "Mooncake" + DifferentiationInterfacePolyesterForwardDiffExt = ["PolyesterForwardDiff", "ForwardDiff", "DiffResults"] + DifferentiationInterfaceReverseDiffExt = ["ReverseDiff", "DiffResults"] + DifferentiationInterfaceSparseArraysExt = "SparseArrays" + DifferentiationInterfaceSparseConnectivityTracerExt = "SparseConnectivityTracer" + DifferentiationInterfaceSparseMatrixColoringsExt = "SparseMatrixColorings" + DifferentiationInterfaceStaticArraysExt = "StaticArrays" + DifferentiationInterfaceSymbolicsExt = "Symbolics" + DifferentiationInterfaceTrackerExt = "Tracker" + DifferentiationInterfaceZygoteExt = ["Zygote", "ForwardDiff"] + + [deps.DifferentiationInterface.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" + Diffractor = "9f5e2b26-1114-432f-b630-d3fe2085c51c" + Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" + EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869" + FastDifferentiation = "eb9bf01b-bf85-4b60-bf87-ee5de06c00be" + FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" + FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" + ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" + GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" + GTPSA = "b27dd330-f138-47c5-815b-40db9dd9b6e8" + Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6" + PolyesterForwardDiff = "98d1487c-24ca-40b6-b7ab-df2af84e126b" + ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" + SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5" + SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35" + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" + Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" + Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" + [[deps.Distances]] deps = ["LinearAlgebra", "Statistics", "StatsAPI"] git-tree-sha1 = "c7e3a542b999843086e2f29dac96a618c105be1d" @@ -380,9 +446,9 @@ version = "1.11.0" [[deps.Distributions]] deps = ["AliasTables", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns"] -git-tree-sha1 = "03aa5d44647eaec98e1920635cdfed5d5560a8b9" +git-tree-sha1 = "3e6d038b77f22791b8e3472b7c633acea1ecac06" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" -version = "0.25.117" +version = "0.25.120" [deps.Distributions.extensions] DistributionsChainRulesCoreExt = "ChainRulesCore" @@ -395,16 +461,15 @@ version = "0.25.117" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [[deps.DocStringExtensions]] -deps = ["LibGit2"] -git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d" +git-tree-sha1 = "7442a5dfe1ebb773c29cc2962a8980f47221d76c" uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.9.3" +version = "0.9.5" [[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 = "d0ea2c044963ed6f37703cead7e29f70cba13d7e" +deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "CodecZlib", "Dates", "DocStringExtensions", "Downloads", "Git", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "TOML", "Test", "Unicode"] +git-tree-sha1 = "e25bc156a7e72f0a9f738815a4426dd01a7e914b" uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -version = "1.8.0" +version = "1.13.0" [[deps.Downloads]] deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] @@ -424,16 +489,21 @@ uuid = "2904ab23-551e-5aed-883f-487f97af5226" version = "0.2.3" [[deps.EnsembleKalmanProcesses]] -deps = ["Convex", "Distributions", "DocStringExtensions", "GaussianRandomFields", "Interpolations", "LinearAlgebra", "MathOptInterface", "Optim", "QuadGK", "Random", "RecipesBase", "SCS", "SparseArrays", "Statistics", "StatsBase", "TOML"] -git-tree-sha1 = "00bb94ff704d7aeed9c72d4a2a05d6abf6cb7946" +deps = ["Convex", "Distributions", "DocStringExtensions", "FFMPEG", "GaussianRandomFields", "Interpolations", "LinearAlgebra", "MathOptInterface", "Optim", "QuadGK", "Random", "RecipesBase", "SCS", "SparseArrays", "Statistics", "StatsBase", "TOML", "TSVD"] +git-tree-sha1 = "dd6b4b258beeef8db8024e14fea95d34305e5bb0" uuid = "aa8a2aa5-91d8-4396-bcef-d4f2ec43552d" -version = "2.0.1" +version = "2.4.0" + +[[deps.EnumX]] +git-tree-sha1 = "bddad79635af6aec424f53ed8aad5d7555dc6f00" +uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56" +version = "1.0.5" [[deps.Expat_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "e51db81749b0777b2147fbe7b783ee79045b8e99" +git-tree-sha1 = "d55dffd9ae73ff72f1c0482454dcf2ec6c6c4a63" uuid = "2e619515-83b5-522b-bb60-26c02a35a201" -version = "2.6.4+3" +version = "2.6.5+0" [[deps.FFMPEG]] deps = ["FFMPEG_jll"] @@ -449,15 +519,15 @@ version = "4.4.4+1" [[deps.FFTW]] deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"] -git-tree-sha1 = "7de7c78d681078f027389e067864a8d53bd7c3c9" +git-tree-sha1 = "797762812ed063b9b94f6cc7742bc8883bb5e69e" uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" -version = "1.8.1" +version = "1.9.0" [[deps.FFTW_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "4d81ed14783ec49ce9f2e168208a12ce1815aa25" +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "6d6219a004b8cf1e0b4dbe27a2860b8e04eba0be" uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a" -version = "3.3.10+3" +version = "3.3.11+0" [[deps.FastGaussQuadrature]] deps = ["LinearAlgebra", "SpecialFunctions", "StaticArrays"] @@ -483,9 +553,9 @@ weakdeps = ["PDMats", "SparseArrays", "Statistics"] [[deps.FiniteDiff]] deps = ["ArrayInterface", "LinearAlgebra", "Setfield"] -git-tree-sha1 = "84e3a47db33be7248daa6274b287507dd6ff84e8" +git-tree-sha1 = "f089ab1f834470c525562030c8cfde4025d5e915" uuid = "6a86dc24-6348-571c-b903-95158fe2bd41" -version = "2.26.2" +version = "2.27.0" [deps.FiniteDiff.extensions] FiniteDiffBandedMatricesExt = "BandedMatrices" @@ -501,9 +571,9 @@ version = "2.26.2" [[deps.Fontconfig_jll]] deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Zlib_jll"] -git-tree-sha1 = "21fac3c77d7b5a9fc03b0ec503aa1a6392c34d2b" +git-tree-sha1 = "301b5d5d731a0654825f1f2e906990f7141a106b" uuid = "a3f928ae-7b40-5064-980b-68af3947d34b" -version = "2.15.0+0" +version = "2.16.0+0" [[deps.Formatting]] deps = ["Logging", "Printf"] @@ -523,15 +593,15 @@ weakdeps = ["StaticArrays"] [[deps.FreeType2_jll]] deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Zlib_jll"] -git-tree-sha1 = "786e968a8d2fb167f2e4880baba62e0e26bd8e4e" +git-tree-sha1 = "2c5512e11c791d1baed2049c5652441b28fc6a31" uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7" -version = "2.13.3+1" +version = "2.13.4+0" [[deps.FriBidi_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "846f7026a9decf3679419122b49f8a1fdb48d2d5" +git-tree-sha1 = "7a214fdac5ed5f59a22c2d9a885a16da1c74bbc7" uuid = "559328eb-81f9-559d-9380-de523a88c83c" -version = "1.0.16+0" +version = "1.0.17+0" [[deps.FunctionWrappers]] git-tree-sha1 = "d62485945ce5ae9c0c48f124a84998d755bae00e" @@ -539,10 +609,10 @@ uuid = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e" version = "1.1.3" [[deps.Functors]] -deps = ["LinearAlgebra"] -git-tree-sha1 = "64d8e93700c7a3f28f717d265382d52fac9fa1c1" +deps = ["Compat", "ConstructionBase", "LinearAlgebra", "Random"] +git-tree-sha1 = "60a0339f28a233601cb74468032b5c302d5067de" uuid = "d9f16b24-f501-4c13-a1f2-28368ffc5196" -version = "0.4.12" +version = "0.5.2" [[deps.Future]] deps = ["Random"] @@ -561,41 +631,41 @@ git-tree-sha1 = "d9c335f2c06424029b2addf9abf602e0feb2f53e" uuid = "e4b2fa32-6e09-5554-b718-106ed5adafe9" version = "2.1.6" -[[deps.Gettext_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"] -git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046" -uuid = "78b55507-aeef-58d4-861c-77aaff3498b1" -version = "0.21.0+0" +[[deps.GettextRuntime_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll"] +git-tree-sha1 = "45288942190db7c5f760f59c04495064eedf9340" +uuid = "b0724c58-0f36-5564-988d-3bb0596ebc4a" +version = "0.22.4+0" [[deps.Git]] -deps = ["Git_jll"] -git-tree-sha1 = "04eff47b1354d702c3a85e8ab23d539bb7d5957e" +deps = ["Git_jll", "JLLWrappers", "OpenSSH_jll"] +git-tree-sha1 = "2230a9cc32394b11a3b3aa807a382e3bbab1198c" uuid = "d7ba0133-e1db-5d97-8f8c-041e4b3a1eb2" -version = "1.3.1" +version = "1.4.0" [[deps.Git_jll]] deps = ["Artifacts", "Expat_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Libiconv_jll", "OpenSSL_jll", "PCRE2_jll", "Zlib_jll"] -git-tree-sha1 = "399f4a308c804b446ae4c91eeafadb2fe2c54ff9" +git-tree-sha1 = "b981ed24de5855f20fce5b8cb767c179f93e4268" uuid = "f8c6e375-362e-5223-8a59-34ff63f689eb" -version = "2.47.1+0" +version = "2.50.0+0" [[deps.Glib_jll]] -deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Zlib_jll"] -git-tree-sha1 = "b0036b392358c80d2d2124746c2bf3d48d457938" +deps = ["Artifacts", "GettextRuntime_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Zlib_jll"] +git-tree-sha1 = "35fbd0cefb04a516104b8e183ce0df11b70a3f1a" uuid = "7746bdde-850d-59dc-9ae8-88ece973131d" -version = "2.82.4+0" +version = "2.84.3+0" [[deps.Graphite2_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "01979f9b37367603e2848ea225918a3b3861b606" +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "8a6dbda1fd736d60cc477d99f2e7a042acfa46e8" uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472" -version = "1.3.14+1" +version = "1.3.15+0" [[deps.HarfBuzz_jll]] deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll"] -git-tree-sha1 = "55c53be97790242c29031e5cd45e8ac296dadda3" +git-tree-sha1 = "f923f9a774fcf3f5cb761bfa43aeadd689714813" uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566" -version = "8.5.0+0" +version = "8.5.1+0" [[deps.HostCPUFeatures]] deps = ["BitTwiddlingConvenienceFunctions", "IfElse", "Libdl", "Static"] @@ -620,9 +690,9 @@ uuid = "22cec73e-a1b8-11e9-2c92-598750a2cf9c" version = "0.3.1" [[deps.InlineStrings]] -git-tree-sha1 = "45521d31238e87ee9f9732561bfee12d4eebd52d" +git-tree-sha1 = "8594fac023c5ce1ef78260f24d1ad18b4327b420" uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48" -version = "1.4.2" +version = "1.4.4" [deps.InlineStrings.extensions] ArrowTypesExt = "ArrowTypes" @@ -656,9 +726,9 @@ version = "0.15.1" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [[deps.IntervalSets]] -git-tree-sha1 = "dba9ddf07f77f60450fe5d2e2beb9854d9a49bd0" +git-tree-sha1 = "5fbb102dcb8b1a858111ae81d56682376130517d" uuid = "8197267c-284f-5f27-9208-e0e47529a953" -version = "0.7.10" +version = "0.7.11" weakdeps = ["Random", "RecipesBase", "Statistics"] [deps.IntervalSets.extensions] @@ -710,9 +780,9 @@ version = "0.21.4" [[deps.JSON3]] deps = ["Dates", "Mmap", "Parsers", "PrecompileTools", "StructTypes", "UUIDs"] -git-tree-sha1 = "1d322381ef7b087548321d3f878cb4c9bd8f8f9b" +git-tree-sha1 = "411eccfe8aba0814ffa0fdf4860913ed09c34975" uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" -version = "1.14.1" +version = "1.14.3" [deps.JSON3.extensions] JSON3ArrowExt = ["ArrowTypes"] @@ -728,15 +798,15 @@ 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" +git-tree-sha1 = "0b8ef8b51580b0d87d0b7a5233bb8ea6d948feb4" uuid = "ec8451be-7e33-11e9-00cf-bbf324bd1392" -version = "0.10.64" +version = "0.10.65" [[deps.LAME_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "170b660facf5df5de098d866564877e119141cbd" +git-tree-sha1 = "059aabebaa7c82ccb853dd4a0ee9d17796f7e1bc" uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d" -version = "3.100.2+0" +version = "3.100.3+0" [[deps.LDLFactorizations]] deps = ["AMD", "LinearAlgebra", "SparseArrays", "Test"] @@ -746,9 +816,9 @@ version = "0.10.1" [[deps.LLVMOpenMP_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "78211fb6cbc872f77cad3fc0b6cf647d923f4929" +git-tree-sha1 = "eb62a3deb62fc6d8822c0c4bef73e4412419c5d8" uuid = "1d63c593-3942-5779-bab2-d838dc0a180e" -version = "18.1.7+0" +version = "18.1.8+0" [[deps.LZO_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -813,22 +883,10 @@ uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" version = "1.11.0" [[deps.Libffi_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "27ecae93dd25ee0909666e6835051dd684cc035e" -uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490" -version = "3.2.2+2" - -[[deps.Libgcrypt_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll"] -git-tree-sha1 = "8be878062e0ffa2c3f67bb58a595375eda5de80b" -uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4" -version = "1.11.0+0" - -[[deps.Libgpg_error_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "df37206100d39f79b3376afb6b9cee4970041c61" -uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8" -version = "1.51.1+0" +git-tree-sha1 = "c8da7e6a91781c41a863611c7e966098d783c57a" +uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490" +version = "3.4.7+0" [[deps.Libiconv_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -838,21 +896,21 @@ version = "1.18.0+0" [[deps.Libmount_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "89211ea35d9df5831fca5d33552c02bd33878419" +git-tree-sha1 = "a31572773ac1b745e0343fe5e2c8ddda7a37e997" uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9" -version = "2.40.3+0" +version = "2.41.0+0" [[deps.Libuuid_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "e888ad02ce716b319e6bdb985d2ef300e7089889" +git-tree-sha1 = "321ccef73a96ba828cd51f2ab5b9f917fa73945a" uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700" -version = "2.40.3+0" +version = "2.41.0+0" [[deps.LineSearches]] deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"] -git-tree-sha1 = "e4c3be53733db1051cc15ecf573b1042b3a712a1" +git-tree-sha1 = "4adee99b7262ad2a1a4bbbc59d993d24e55ea96f" uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" -version = "7.3.0" +version = "7.4.0" [[deps.LinearAlgebra]] deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] @@ -893,9 +951,9 @@ version = "1.1.0" [[deps.LoopVectorization]] deps = ["ArrayInterface", "CPUSummary", "CloseOpenIntervals", "DocStringExtensions", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "PrecompileTools", "SIMDTypes", "SLEEFPirates", "Static", "StaticArrayInterface", "ThreadingUtilities", "UnPack", "VectorizationBase"] -git-tree-sha1 = "8084c25a250e00ae427a379a5b607e7aed96a2dd" +git-tree-sha1 = "e5afce7eaf5b5ca0d444bcb4dc4fd78c54cbbac0" uuid = "bdcacae8-1622-11e9-2a5c-532679323890" -version = "0.12.171" +version = "0.12.172" weakdeps = ["ChainRulesCore", "ForwardDiff", "SpecialFunctions"] [deps.LoopVectorization.extensions] @@ -916,20 +974,20 @@ version = "0.2.1" [[deps.MKL_jll]] deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "oneTBB_jll"] -git-tree-sha1 = "ed4097130e3dd3721814b5f277da72f48905e80c" +git-tree-sha1 = "5de60bc6cb3899cd318d80d627560fae2e2d99ae" uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7" -version = "2025.0.1+0" +version = "2025.0.1+1" [[deps.MLJModelInterface]] -deps = ["Random", "ScientificTypesBase", "StatisticalTraits"] -git-tree-sha1 = "ceaff6618408d0e412619321ae43b33b40c1a733" +deps = ["REPL", "Random", "ScientificTypesBase", "StatisticalTraits"] +git-tree-sha1 = "66626f80d5807921045d539b4f7153b1d47c5f8a" uuid = "e80e1ace-859a-464e-9ed9-23947d8ae3ea" -version = "1.11.0" +version = "1.11.1" [[deps.MacroTools]] -git-tree-sha1 = "72aebe0b5051e5143a079a4685a46da330a40472" +git-tree-sha1 = "1e0228a030642014fe5cfe68c2c0a818f9e3f522" uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" -version = "0.5.15" +version = "0.5.16" [[deps.ManualMemory]] git-tree-sha1 = "bcaef4fc7a0cfe2cba636d84cda54b5e4e4ca3cd" @@ -948,10 +1006,10 @@ uuid = "d0879d2d-cac2-40c8-9cee-1863dc0c7391" version = "0.1.2" [[deps.MathOptInterface]] -deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "DataStructures", "ForwardDiff", "JSON3", "LinearAlgebra", "MutableArithmetics", "NaNMath", "OrderedCollections", "PrecompileTools", "Printf", "SparseArrays", "SpecialFunctions", "Test", "Unicode"] -git-tree-sha1 = "f346e3b50c8bb62b7a5c4961ebdbaef65d1ce1e2" +deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "DataStructures", "ForwardDiff", "JSON3", "LinearAlgebra", "MutableArithmetics", "NaNMath", "OrderedCollections", "PrecompileTools", "Printf", "SparseArrays", "SpecialFunctions", "Test"] +git-tree-sha1 = "2f2c18c6acab9042557bdb0af8c3a14dd7b64413" uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" -version = "1.35.1" +version = "1.41.0" [[deps.MbedTLS_jll]] deps = ["Artifacts", "Libdl"] @@ -980,21 +1038,21 @@ version = "2023.12.12" [[deps.MutableArithmetics]] deps = ["LinearAlgebra", "SparseArrays", "Test"] -git-tree-sha1 = "43122df26d27424b23577d59e2d8020f28386516" +git-tree-sha1 = "491bdcdc943fcbc4c005900d7463c9f216aabf4c" uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" -version = "1.6.2" +version = "1.6.4" [[deps.NLSolversBase]] -deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"] -git-tree-sha1 = "a0b464d183da839699f4c79e7606d9d186ec172c" +deps = ["ADTypes", "DifferentiationInterface", "Distributed", "FiniteDiff", "ForwardDiff"] +git-tree-sha1 = "25a6638571a902ecfb1ae2a18fc1575f86b1d4df" uuid = "d41bc354-129a-5804-8e4c-c37616107c6c" -version = "7.8.3" +version = "7.10.0" [[deps.NaNMath]] deps = ["OpenLibm_jll"] -git-tree-sha1 = "fe891aea7ccd23897520db7f16931212454e277e" +git-tree-sha1 = "9b8215b1ee9e78a293f99797cd31375471b2bcae" uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" -version = "1.1.1" +version = "1.1.3" [[deps.NaturalSort]] git-tree-sha1 = "eda490d06b9f7c00752ee81cfa451efe55521e21" @@ -1006,19 +1064,19 @@ uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" version = "1.2.0" [[deps.OffsetArrays]] -git-tree-sha1 = "5e1897147d1ff8d98883cda2be2187dcf57d8f0c" +git-tree-sha1 = "117432e406b5c023f665fa73dc26e79ec3630151" uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" -version = "1.15.0" +version = "1.17.0" weakdeps = ["Adapt"] [deps.OffsetArrays.extensions] OffsetArraysAdaptExt = "Adapt" [[deps.Ogg_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "887579a3eb005446d514ab7aeac5d1d027658b8f" +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "b6aa4566bb7ae78498a5e68943863fa8b5231b59" uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051" -version = "1.3.5+1" +version = "1.3.6+0" [[deps.OpenBLAS32_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"] @@ -1034,13 +1092,19 @@ version = "0.3.27+1" [[deps.OpenLibm_jll]] deps = ["Artifacts", "Libdl"] uuid = "05823500-19ac-5b8b-9628-191a04bc5112" -version = "0.8.1+2" +version = "0.8.5+0" + +[[deps.OpenSSH_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "OpenSSL_jll", "Zlib_jll"] +git-tree-sha1 = "cb7acd5d10aff809b4d0191dfe1956c2edf35800" +uuid = "9bd350c2-7e96-507f-8002-3f2e150b4e1b" +version = "10.0.1+0" [[deps.OpenSSL_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "7493f61f55a6cce7325f197443aa80d32554ba10" +git-tree-sha1 = "87510f7292a2b21aeff97912b0898f9553cc5c2c" uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" -version = "3.0.15+3" +version = "3.5.1+0" [[deps.OpenSpecFun_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"] @@ -1049,10 +1113,10 @@ uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" version = "0.5.6+0" [[deps.Optim]] -deps = ["Compat", "FillArrays", "ForwardDiff", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"] -git-tree-sha1 = "ab7edad78cdef22099f43c54ef77ac63c2c9cc64" +deps = ["Compat", "EnumX", "FillArrays", "ForwardDiff", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"] +git-tree-sha1 = "61942645c38dd2b5b78e2082c9b51ab315315d10" uuid = "429524aa-4258-5aef-a3af-852621145aeb" -version = "1.10.0" +version = "1.13.2" weakdeps = ["MathOptInterface"] [deps.Optim.extensions] @@ -1060,14 +1124,14 @@ weakdeps = ["MathOptInterface"] [[deps.Opus_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "6703a85cb3781bd5909d48730a67205f3f31a575" +git-tree-sha1 = "c392fc5dd032381919e3b22dd32d6443760ce7ea" uuid = "91d4177d-7536-5919-b921-800302f37372" -version = "1.3.3+0" +version = "1.5.2+0" [[deps.OrderedCollections]] -git-tree-sha1 = "12f1439c4f986bb868acda6ea33ebc78e19b95ad" +git-tree-sha1 = "05868e21324cede2207c6f0f466b4bfef6d5e7ee" uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.7.0" +version = "1.8.1" [[deps.PCRE2_jll]] deps = ["Artifacts", "Libdl"] @@ -1076,9 +1140,9 @@ version = "10.42.0+1" [[deps.PDMats]] deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"] -git-tree-sha1 = "966b85253e959ea89c53a9abebbf2e964fbf593b" +git-tree-sha1 = "f07c06228a1c670ae4c87d1276b92c7c597fdda0" uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" -version = "0.11.32" +version = "0.11.35" [[deps.Parameters]] deps = ["OrderedCollections", "UnPack"] @@ -1088,15 +1152,15 @@ version = "0.12.3" [[deps.Parsers]] deps = ["Dates", "PrecompileTools", "UUIDs"] -git-tree-sha1 = "8489905bcdbcfac64d1daa51ca07c0d8f0283821" +git-tree-sha1 = "7d2f8f21da5db6a806faf7b9b292296da42b2810" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "2.8.1" +version = "2.8.3" [[deps.Pixman_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl"] -git-tree-sha1 = "35621f10a7531bc8fa58f74610b1bfb70a3cfc6b" +git-tree-sha1 = "db76b1ecd5e9715f3d043cec13b2ec93ce015d53" uuid = "30392449-352a-5448-841d-b1acce4e97dc" -version = "0.43.4+0" +version = "0.44.2+0" [[deps.Pkg]] deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "Random", "SHA", "TOML", "Tar", "UUIDs", "p7zip_jll"] @@ -1160,15 +1224,15 @@ version = "1.5.1" [[deps.ProgressLogging]] deps = ["Logging", "SHA", "UUIDs"] -git-tree-sha1 = "80d919dee55b9c50e8d9e2da5eeafff3fe58b539" +git-tree-sha1 = "d95ed0324b0799843ac6f7a6a85e65fe4e5173f0" uuid = "33c8b6b6-d38a-422a-b730-caa89a2f386c" -version = "0.1.4" +version = "0.1.5" [[deps.ProgressMeter]] deps = ["Distributed", "Printf"] -git-tree-sha1 = "8f6bc219586aef8baf0ff9a5fe16ee9c70cb65e4" +git-tree-sha1 = "13c5103482a8ed1536a54c08d0e742ae3dca2d42" uuid = "92933f4c-e287-5a05-a399-4b506db050ca" -version = "1.10.2" +version = "1.10.4" [[deps.PtrArrays]] git-tree-sha1 = "1d36ef11a9aaf1e8b74dacc6a731dd1de8fd493d" @@ -1183,9 +1247,9 @@ version = "1.96.4" [[deps.QuadGK]] deps = ["DataStructures", "LinearAlgebra"] -git-tree-sha1 = "cda3b045cf9ef07a08ad46731f5a3165e56cf3da" +git-tree-sha1 = "9da16da70037ba9d701192e27befedefb91ec284" uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" -version = "2.11.1" +version = "2.11.2" [deps.QuadGK.extensions] QuadGKEnzymeExt = "Enzyme" @@ -1245,15 +1309,15 @@ version = "0.1.0" [[deps.Requires]] deps = ["UUIDs"] -git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7" +git-tree-sha1 = "62389eeff14780bfe55195b7204c0d8738436d64" uuid = "ae029012-a4dd-5104-9daa-d747884805df" -version = "1.3.0" +version = "1.3.1" [[deps.ReverseDiff]] deps = ["ChainRulesCore", "DiffResults", "DiffRules", "ForwardDiff", "FunctionWrappers", "LinearAlgebra", "LogExpFunctions", "MacroTools", "NaNMath", "Random", "SpecialFunctions", "StaticArrays", "Statistics"] -git-tree-sha1 = "cc6cd622481ea366bb9067859446a8b01d92b468" +git-tree-sha1 = "3ab8eee3620451b09f0272c271875b4bc02146d9" uuid = "37e2e3b7-166d-5795-8a7a-e32c996b4267" -version = "1.15.3" +version = "1.16.1" [[deps.Rmath]] deps = ["Random", "Rmath_jll"] @@ -1268,10 +1332,10 @@ uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f" version = "0.4.3+0" [[deps.SCS]] -deps = ["MathOptInterface", "Requires", "SCS_jll", "SparseArrays"] -git-tree-sha1 = "aa3fcff53da363b4ba4b54d4ac4c9186ab00d703" +deps = ["MathOptInterface", "PrecompileTools", "SCS_jll", "SparseArrays"] +git-tree-sha1 = "4aed85dec0209d638c241c34160999eaaf07965a" uuid = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13" -version = "2.0.2" +version = "2.1.0" [deps.SCS.extensions] SCSSCS_GPU_jllExt = ["SCS_GPU_jll"] @@ -1331,9 +1395,9 @@ version = "1.11.0" [[deps.Setfield]] deps = ["ConstructionBase", "Future", "MacroTools", "StaticArraysCore"] -git-tree-sha1 = "e2cc6d8c88613c05e1defb55170bf5ff211fbeac" +git-tree-sha1 = "c5391c6ace3bc430ca630251d02ea9687169ca68" uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46" -version = "1.1.1" +version = "1.1.2" [[deps.SharedArrays]] deps = ["Distributed", "Mmap", "Random", "Serialization"] @@ -1357,9 +1421,9 @@ version = "1.11.0" [[deps.SpecialFunctions]] deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] -git-tree-sha1 = "64cca0c26b4f31ba18f13f6c12af7c85f478cfde" +git-tree-sha1 = "41852b8679f78c8d8961eeadc8f62cef861a52e3" uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "2.5.0" +version = "2.5.1" weakdeps = ["ChainRulesCore"] [deps.SpecialFunctions.extensions] @@ -1373,15 +1437,15 @@ version = "0.1.15" [[deps.StableRNGs]] deps = ["Random"] -git-tree-sha1 = "83e6cce8324d49dfaf9ef059227f91ed4441a8e5" +git-tree-sha1 = "95af145932c2ed859b63329952ce8d633719f091" uuid = "860ef19b-820b-49d6-a774-d7a799459cd3" -version = "1.0.2" +version = "1.0.3" [[deps.Static]] deps = ["CommonWorldInvalidations", "IfElse", "PrecompileTools"] -git-tree-sha1 = "87d51a3ee9a4b0d2fe054bdd3fc2436258db2603" +git-tree-sha1 = "f737d444cb0ad07e61b3c1bef8eb91203c321eff" uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" -version = "1.1.1" +version = "1.2.0" [[deps.StaticArrayInterface]] deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "Static"] @@ -1396,9 +1460,9 @@ weakdeps = ["OffsetArrays", "StaticArrays"] [[deps.StaticArrays]] deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"] -git-tree-sha1 = "47091a0340a675c738b1304b58161f3b0839d454" +git-tree-sha1 = "0feb6b9031bd5c51f9072393eb5ab3efd31bf9e4" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.9.10" +version = "1.9.13" weakdeps = ["ChainRulesCore", "Statistics"] [deps.StaticArrays.extensions] @@ -1428,9 +1492,9 @@ weakdeps = ["SparseArrays"] [[deps.StatsAPI]] deps = ["LinearAlgebra"] -git-tree-sha1 = "1ff449ad350c9c4cbc756624d6f8a8c3ef56d3ed" +git-tree-sha1 = "9d72a13a3f4dd3795a195ac5a44d7d6ff5f552ff" uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" -version = "1.7.0" +version = "1.7.1" [[deps.StatsBase]] deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"] @@ -1446,9 +1510,9 @@ version = "0.9.18" [[deps.StringManipulation]] deps = ["PrecompileTools"] -git-tree-sha1 = "a6b1675a536c5ad1a60e5a5153e1fee12eb146e3" +git-tree-sha1 = "725421ae8e530ec29bcbdddbe91ff8053421d023" uuid = "892a3eda-7b42-436c-8928-eab12a02cf0e" -version = "0.4.0" +version = "0.4.1" [[deps.StructTypes]] deps = ["Dates", "UUIDs"] @@ -1474,6 +1538,12 @@ deps = ["Dates"] uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" version = "1.0.3" +[[deps.TSVD]] +deps = ["Adapt", "LinearAlgebra"] +git-tree-sha1 = "c39caef6bae501e5607a6caf68dd9ac6e8addbcb" +uuid = "9449cd9e-2762-5aa3-a617-5413e99d722e" +version = "0.4.4" + [[deps.TableTraits]] deps = ["IteratorInterfaceExtensions"] git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39" @@ -1482,9 +1552,9 @@ version = "1.0.1" [[deps.Tables]] deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "OrderedCollections", "TableTraits"] -git-tree-sha1 = "598cd7c1f68d1e205689b1c2fe65a9f85846f297" +git-tree-sha1 = "f2c1efbc8f3a609aadf318094f8fc5204bdaf344" uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" -version = "1.12.0" +version = "1.12.1" [[deps.Tar]] deps = ["ArgTools", "SHA"] @@ -1510,9 +1580,9 @@ version = "1.11.0" [[deps.ThreadingUtilities]] deps = ["ManualMemory"] -git-tree-sha1 = "eda08f7e9818eb53661b3deb74e3159460dfbc27" +git-tree-sha1 = "d969183d3d244b6c33796b5ed01ab97328f2db85" uuid = "8290d209-cae3-49c0-8002-c8c24d57dab5" -version = "0.5.2" +version = "0.5.5" [[deps.TranscodingStreams]] git-tree-sha1 = "0c45878dcfdcfa8480052b6ab162cdd138781742" @@ -1541,9 +1611,9 @@ version = "0.4.80" [[deps.Tullio]] deps = ["DiffRules", "LinearAlgebra", "Requires"] -git-tree-sha1 = "6d476962ba4e435d7f4101a403b1d3d72afe72f3" +git-tree-sha1 = "972698b132b9df8791ae74aa547268e977b55f68" uuid = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc" -version = "0.3.7" +version = "0.3.8" [deps.Tullio.extensions] TullioCUDAExt = "CUDA" @@ -1588,65 +1658,47 @@ git-tree-sha1 = "c1a7aa6219628fcd757dede0ca95e245c5cd9511" uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6" version = "1.0.0" -[[deps.XML2_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"] -git-tree-sha1 = "a2fccc6559132927d4c5dc183e3e01048c6dcbd6" -uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a" -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 = "7d1671acbe47ac88e981868a078bd6b4e27c5191" -uuid = "aed1982a-8fda-507f-9586-7b0439959a61" -version = "1.1.42+0" - [[deps.Xorg_libX11_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libxcb_jll", "Xorg_xtrans_jll"] -git-tree-sha1 = "9dafcee1d24c4f024e7edc92603cedba72118283" +git-tree-sha1 = "b5899b25d17bf1889d25906fb9deed5da0c15b3b" uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc" -version = "1.8.6+3" +version = "1.8.12+0" [[deps.Xorg_libXau_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "e9216fdcd8514b7072b43653874fd688e4c6c003" +git-tree-sha1 = "aa1261ebbac3ccc8d16558ae6799524c450ed16b" uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec" -version = "1.0.12+0" +version = "1.0.13+0" [[deps.Xorg_libXdmcp_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "89799ae67c17caa5b3b5a19b8469eeee474377db" +git-tree-sha1 = "52858d64353db33a56e13c341d7bf44cd0d7b309" uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05" -version = "1.1.5+0" +version = "1.1.6+0" [[deps.Xorg_libXext_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll"] -git-tree-sha1 = "d7155fea91a4123ef59f42c4afb5ab3b4ca95058" +git-tree-sha1 = "a4c0ee07ad36bf8bbce1c3bb52d21fb1e0b987fb" uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3" -version = "1.3.6+3" +version = "1.3.7+0" [[deps.Xorg_libXrender_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll"] -git-tree-sha1 = "a490c6212a0e90d2d55111ac956f7c4fa9c277a6" +git-tree-sha1 = "7ed9347888fac59a618302ee38216dd0379c480d" uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa" -version = "0.9.11+1" - -[[deps.Xorg_libpthread_stubs_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "c57201109a9e4c0585b208bb408bc41d205ac4e9" -uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74" -version = "0.1.2+0" +version = "0.9.12+0" [[deps.Xorg_libxcb_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"] -git-tree-sha1 = "1a74296303b6524a0472a8cb12d3d87a78eb3612" +deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libXau_jll", "Xorg_libXdmcp_jll"] +git-tree-sha1 = "bfcaf7ec088eaba362093393fe11aa141fa15422" uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b" -version = "1.17.0+3" +version = "1.17.1+0" [[deps.Xorg_xtrans_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "6dba04dbfb72ae3ebe5418ba33d087ba8aa8cb00" +git-tree-sha1 = "a63799ff68005991f9d9491b6e95bd3478d783cb" uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10" -version = "1.5.1+0" +version = "1.6.0+0" [[deps.Zlib_jll]] deps = ["Libdl"] @@ -1678,21 +1730,21 @@ version = "5.11.0+0" [[deps.libfdk_aac_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "8a22cf860a7d27e4f3498a0fe0811a7957badb38" +git-tree-sha1 = "646634dd19587a56ee2f1199563ec056c5f228df" uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280" -version = "2.0.3+0" +version = "2.0.4+0" [[deps.libpng_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Zlib_jll"] -git-tree-sha1 = "d7b5bbf1efbafb5eca466700949625e07533aff2" +git-tree-sha1 = "07b6a107d926093898e82b3b1db657ebe33134ec" uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f" -version = "1.6.45+1" +version = "1.6.50+0" [[deps.libvorbis_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"] -git-tree-sha1 = "490376214c4721cdaca654041f635213c6165cb3" +deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll"] +git-tree-sha1 = "11e1772e7f3cc987e9d3de991dd4f6b2602663a5" uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a" -version = "1.3.7+2" +version = "1.3.8+0" [[deps.nghttp2_jll]] deps = ["Artifacts", "Libdl"] @@ -1701,9 +1753,9 @@ version = "1.59.0+0" [[deps.oneTBB_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "7d0ea0f4895ef2f5cb83645fa689e52cb55cf493" +git-tree-sha1 = "d5a767a3bb77135a99e433afe0eb14cd7f6914c3" uuid = "1317d2d5-d96f-522e-a858-c73665f53c3e" -version = "2021.12.0+0" +version = "2022.0.0+0" [[deps.p7zip_jll]] deps = ["Artifacts", "Libdl"] diff --git a/examples/Cloudy/Cloudy_calibrate.jl b/examples/Cloudy/Cloudy_calibrate.jl index 1b4fd34b5..f3fb84db0 100644 --- a/examples/Cloudy/Cloudy_calibrate.jl +++ b/examples/Cloudy/Cloudy_calibrate.jl @@ -94,8 +94,8 @@ prior_N0 = constrained_gaussian(param_names[1], 400, 300, 0.4 * N0_true, Inf) prior_θ = constrained_gaussian(param_names[2], 1.0, 5.0, 1e-1, Inf) prior_k = constrained_gaussian(param_names[3], 0.2, 1.0, 1e-4, Inf) priors = combine_distributions([prior_N0, prior_θ, prior_k]) -# Plot the priors -p = plot(priors, constrained = false) +# Plots.Plot the priors +p = Plots.plot(priors, constrained = false) savefig(p, output_directory * "cloudy_priors.png") ### @@ -160,7 +160,7 @@ dummy = ones(n_params) dist_type = ParticleDistributions.GammaPrimitiveParticleDistribution(dummy...) model_settings = DynamicalModel.ModelSettings(kernel, dist_type, moments, tspan) # EKI iterations -n_iter = N_iter +n_iter = [N_iter] for n in 1:N_iter # Return transformed parameters in physical/constrained space ϕ_n = get_ϕ_final(priors, ekiobj) @@ -169,12 +169,12 @@ for n in 1:N_iter G_ens = hcat(G_n...) # reformat terminate = EnsembleKalmanProcesses.update_ensemble!(ekiobj, G_ens) if !isnothing(terminate) - n_iter = n - 1 + n_iter[1] = n - 1 break end end - +n_iter = n_iter[1] # EKI results: Has the ensemble collapsed toward the truth? θ_true = transform_constrained_to_unconstrained(priors, ϕ_true) @@ -211,8 +211,14 @@ u_init = get_u_prior(ekiobj) anim_eki_unconst_cloudy = @animate for i in 1:(n_iter - 1) u_i = get_u(ekiobj, i) - p1 = plot(u_i[1, :], u_i[2, :], seriestype = :scatter, xlims = extrema(u_init[1, :]), ylims = extrema(u_init[2, :])) - plot!( + p1 = Plots.plot( + u_i[1, :], + u_i[2, :], + seriestype = :scatter, + xlims = extrema(u_init[1, :]), + ylims = extrema(u_init[2, :]), + ) + Plots.plot!( p1, [θ_true[1]], xaxis = "u1", @@ -224,10 +230,16 @@ anim_eki_unconst_cloudy = @animate for i in 1:(n_iter - 1) margin = 5mm, title = "EKI iteration = " * string(i), ) - plot!(p1, [θ_true[2]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = "optimum") - - p2 = plot(u_i[2, :], u_i[3, :], seriestype = :scatter, xlims = extrema(u_init[2, :]), ylims = extrema(u_init[3, :])) - plot!( + Plots.plot!(p1, [θ_true[2]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = "optimum") + + p2 = Plots.plot( + u_i[2, :], + u_i[3, :], + seriestype = :scatter, + xlims = extrema(u_init[2, :]), + ylims = extrema(u_init[3, :]), + ) + Plots.plot!( p2, [θ_true[2]], xaxis = "u2", @@ -240,10 +252,16 @@ anim_eki_unconst_cloudy = @animate for i in 1:(n_iter - 1) title = "EKI iteration = " * string(i), ) - plot!(p2, [θ_true[3]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = "optimum") + Plots.plot!(p2, [θ_true[3]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = "optimum") - p3 = plot(u_i[3, :], u_i[1, :], seriestype = :scatter, xlims = extrema(u_init[3, :]), ylims = extrema(u_init[1, :])) - plot!( + p3 = Plots.plot( + u_i[3, :], + u_i[1, :], + seriestype = :scatter, + xlims = extrema(u_init[3, :]), + ylims = extrema(u_init[1, :]), + ) + Plots.plot!( p3, [θ_true[3]], xaxis = "u3", @@ -256,9 +274,9 @@ anim_eki_unconst_cloudy = @animate for i in 1:(n_iter - 1) title = "EKI iteration = " * string(i), ) - plot!(p3, [θ_true[1]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = "optimum") + Plots.plot!(p3, [θ_true[1]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = "optimum") - p = plot(p1, p2, p3, layout = (1, 3)) + p = Plots.plot(p1, p2, p3, layout = (1, 3)) end gif(anim_eki_unconst_cloudy, joinpath(output_directory, "cloudy_eki_unconstr.gif"), fps = 1) # hide @@ -268,8 +286,14 @@ gif(anim_eki_unconst_cloudy, joinpath(output_directory, "cloudy_eki_unconstr.gif anim_eki_cloudy = @animate for i in 1:(n_iter - 1) ϕ_i = get_ϕ(priors, ekiobj, i) - p1 = plot(ϕ_i[1, :], ϕ_i[2, :], seriestype = :scatter, xlims = extrema(ϕ_init[1, :]), ylims = extrema(ϕ_init[2, :])) - plot!( + p1 = Plots.plot( + ϕ_i[1, :], + ϕ_i[2, :], + seriestype = :scatter, + xlims = extrema(ϕ_init[1, :]), + ylims = extrema(ϕ_init[2, :]), + ) + Plots.plot!( p1, [ϕ_true[1]], xaxis = "ϕ1", @@ -282,11 +306,17 @@ anim_eki_cloudy = @animate for i in 1:(n_iter - 1) title = "EKI iteration = " * string(i), ) - plot!(p1, [ϕ_true[2]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = "optimum") + Plots.plot!(p1, [ϕ_true[2]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = "optimum") - p2 = plot(ϕ_i[2, :], ϕ_i[3, :], seriestype = :scatter, xlims = extrema(ϕ_init[2, :]), ylims = extrema(ϕ_init[3, :])) + p2 = Plots.plot( + ϕ_i[2, :], + ϕ_i[3, :], + seriestype = :scatter, + xlims = extrema(ϕ_init[2, :]), + ylims = extrema(ϕ_init[3, :]), + ) - plot!( + Plots.plot!( p2, [ϕ_true[2]], xaxis = "ϕ2", @@ -299,10 +329,16 @@ anim_eki_cloudy = @animate for i in 1:(n_iter - 1) title = "EKI iteration = " * string(i), ) - plot!(p2, [ϕ_true[3]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = "optimum") + Plots.plot!(p2, [ϕ_true[3]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = "optimum") - p3 = plot(ϕ_i[3, :], ϕ_i[1, :], seriestype = :scatter, xlims = extrema(ϕ_init[3, :]), ylims = extrema(ϕ_init[1, :])) - plot!( + p3 = Plots.plot( + ϕ_i[3, :], + ϕ_i[1, :], + seriestype = :scatter, + xlims = extrema(ϕ_init[3, :]), + ylims = extrema(ϕ_init[1, :]), + ) + Plots.plot!( p3, [ϕ_true[3]], xaxis = "ϕ3", @@ -314,9 +350,9 @@ anim_eki_cloudy = @animate for i in 1:(n_iter - 1) label = false, title = "EKI iteration = " * string(i), ) - plot!(p3, [ϕ_true[1]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = "optimum") + Plots.plot!(p3, [ϕ_true[1]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = "optimum") - p = plot(p1, p2, p3, layout = (1, 3)) + p = Plots.plot(p1, p2, p3, layout = (1, 3)) end diff --git a/examples/Cloudy/Cloudy_emulate_sample.jl b/examples/Cloudy/Cloudy_emulate_sample.jl index c6b04232e..4c1a2b514 100644 --- a/examples/Cloudy/Cloudy_emulate_sample.jl +++ b/examples/Cloudy/Cloudy_emulate_sample.jl @@ -17,6 +17,8 @@ using CalibrateEmulateSample.Emulators using CalibrateEmulateSample.MarkovChainMonteCarlo using CalibrateEmulateSample.Utilities using CalibrateEmulateSample.EnsembleKalmanProcesses +using CalibrateEmulateSample.ParameterDistributions +using CalibrateEmulateSample.DataContainers ################################################################################ # # @@ -90,8 +92,8 @@ function main() cases = [ "rf-scalar", - "gp-gpjl", # Veeeery slow predictions "rf-nonsep", + "gp-gpjl", # Veeeery slow predictions ] # Specify cases to run (e.g., case_mask = [2] only runs the second case) diff --git a/examples/Emulator/G-function/emulate.jl b/examples/Emulator/G-function/emulate.jl index 91ad687fe..6c4afd02e 100644 --- a/examples/Emulator/G-function/emulate.jl +++ b/examples/Emulator/G-function/emulate.jl @@ -40,8 +40,8 @@ function main() rng = MersenneTwister(seed) - n_repeats = 5# repeat exp with same data. - n_dimensions = 10 + n_repeats = 10# repeat exp with same data. + n_dimensions = 20 # To create the sampling n_data_gen = 800 @@ -81,7 +81,8 @@ function main() input[:, i] = samples[ind[i], :] output[i] = y[ind[i]] + noise[i] end - encoder_schedule = (decorrelate_structure_mat(), "out") + # encoder_schedule = (decorrelate_structure_mat(), "out") + encoder_schedule = (minmax_scale(), "in_and_out") iopairs = PairedDataContainer(input, output) # analytic sobol indices taken from @@ -99,12 +100,11 @@ function main() overrides = Dict( "verbose" => true, "scheduler" => DataMisfitController(terminate_at = 1e2), - "n_features_opt" => 120, + "n_features_opt" => 100,#120, "n_iteration" => 10, - "cov_sample_multiplier" => 3.0, - #"localization" => SEC(0.1),#,Doesnt help much tbh + "cov_sample_multiplier" => 0.1, "n_ensemble" => 100, #40*n_dimensions, - "n_cross_val_sets" => 4, + "n_cross_val_sets" => 2, ) if case == "Prior" # don't do anything @@ -126,7 +126,8 @@ function main() elseif case ∈ ["RF-scalar", "Prior"] rank = n_dimensions #<= 10 ? n_dimensions : 10 - kernel_structure = SeparableKernel(LowRankFactor(rank, nugget), OneDimFactor()) + # kernel_structure = SeparableKernel(LowRankFactor(rank, nugget), OneDimFactor()) + kernel_structure = SeparableKernel(DiagonalFactor(nugget), OneDimFactor()) n_features = n_dimensions <= 10 ? n_dimensions * 100 : 1000 if (n_features / n_train_pts > 0.9) && (n_features / n_train_pts < 1.1) @warn "The number of features similar to the number of training points, poor performance expected, change one or other of these" diff --git a/examples/Emulator/G-function/plot_result.jl b/examples/Emulator/G-function/plot_result.jl index 03455b87a..600d73513 100644 --- a/examples/Emulator/G-function/plot_result.jl +++ b/examples/Emulator/G-function/plot_result.jl @@ -7,8 +7,8 @@ function main() output_directory = "output" cases = ["Prior", "GP", "RF-scalar"] case = cases[3] - n_dimensions = 3 - filename = joinpath(output_directory, "Gfunction_$(case)_$(n_dimensions).jld2") + n_dimensions = 20 + filename = joinpath(output_directory, "GFunction_$(case)_$(n_dimensions).jld2") legend = true ( diff --git a/examples/Emulator/Ishigami/emulate.jl b/examples/Emulator/Ishigami/emulate.jl index 105304a48..6136590ac 100644 --- a/examples/Emulator/Ishigami/emulate.jl +++ b/examples/Emulator/Ishigami/emulate.jl @@ -37,7 +37,7 @@ function main() rng = MersenneTwister(seed) - n_repeats = 30 # repeat exp with same data. + n_repeats = 10 # repeat exp with same data. # To create the sampling n_data_gen = 2000 @@ -82,14 +82,15 @@ function main() end iopairs = PairedDataContainer(input, output) cases = ["Prior", "GP", "RF-scalar"] - case = cases[2] + case = cases[3] encoder_schedule = (decorrelate_structure_mat(), "out") - nugget = Float64(1e-4) + nugget = Float64(1e-8) overrides = Dict( + "verbose" => true, "scheduler" => DataMisfitController(terminate_at = 1e4), - "n_features_opt" => 150, + "n_features_opt" => 100, "n_ensemble" => 50, - "n_iteration" => 20, + "n_iteration" => 10, ) if case == "Prior" diff --git a/examples/Emulator/Ishigami/plot_result.jl b/examples/Emulator/Ishigami/plot_result.jl index 18345977e..e36df1cd0 100644 --- a/examples/Emulator/Ishigami/plot_result.jl +++ b/examples/Emulator/Ishigami/plot_result.jl @@ -9,7 +9,7 @@ function main() output_directory = "output" cases = ["Prior", "GP", "RF-scalar"] - case = cases[2] + case = cases[3] filename = joinpath(output_directory, "results_$case.jld2") (sobol_pts, train_idx, mlt_pred_y, mlt_sobol, analytic_sobol, true_y, noise_sample, noise_cov, estimated_sobol) = diff --git a/examples/Emulator/L63/emulate.jl b/examples/Emulator/L63/emulate.jl index 3702c58ec..026e3c51a 100644 --- a/examples/Emulator/L63/emulate.jl +++ b/examples/Emulator/L63/emulate.jl @@ -26,7 +26,7 @@ function main() # rng rng = MersenneTwister(1232435) - n_repeats = 1 # repeat exp with same data. + n_repeats = 5 # repeat exp with same data. println("run experiment $n_repeats times") #for later plots @@ -96,20 +96,11 @@ function main() # Emulate - cases = [ - "GP", - "RF-prior", - "RF-scalar", - "RF-scalar-diagin", - "RF-svd-nonsep", - "RF-nosvd-nonsep", - "RF-nosvd-sep", - "RF-svd-sep", - ] - - case = cases[3] # 7 - - nugget = Float64(1e-8) + cases = ["GP", "RF-prior", "RF-lr-scalar", "rf-diag-scalar", "RF-lr-lr", "RF-nonsep"] + + case = cases[4] # 7 + + nugget = Float64(1e-4) u_test = [] u_hist = [] train_err = [] @@ -148,10 +139,10 @@ function main() optimizer_options = rf_optimizer_overrides, ) - elseif case ∈ ["RF-scalar", "RF-scalar-diagin"] + elseif case ∈ ["RF-lr-scalar", "rf-diag-scalar"] n_features = 10 * Int(floor(sqrt(3 * n_tp))) kernel_structure = - case == "RF-scalar-diagin" ? SeparableKernel(DiagonalFactor(nugget), OneDimFactor()) : + case == "rf-diag-scalar" ? SeparableKernel(DiagonalFactor(nugget), OneDimFactor()) : SeparableKernel(LowRankFactor(3, nugget), OneDimFactor()) mlt = ScalarRandomFeatureInterface( n_features, @@ -160,7 +151,7 @@ function main() kernel_structure = kernel_structure, optimizer_options = rf_optimizer_overrides, ) - elseif case ∈ ["RF-svd-nonsep"] + elseif case ∈ ["RF-nonsep"] kernel_structure = NonseparableKernel(LowRankFactor(4, nugget)) n_features = 500 @@ -172,18 +163,7 @@ function main() kernel_structure = kernel_structure, optimizer_options = rf_optimizer_overrides, ) - elseif case ∈ ["RF-nosvd-nonsep"] - kernel_structure = NonseparableKernel(LowRankFactor(4, nugget)) - n_features = 500 - mlt = VectorRandomFeatureInterface( - n_features, - 3, - 3, - rng = rng, - kernel_structure = kernel_structure, - optimizer_options = rf_optimizer_overrides, - ) - elseif case ∈ ["RF-nosvd-sep", "RF-svd-sep"] + elseif case ∈ ["RF-lr-lr"] kernel_structure = SeparableKernel(LowRankFactor(3, nugget), LowRankFactor(1, nugget)) n_features = 500 mlt = VectorRandomFeatureInterface( @@ -210,7 +190,7 @@ function main() end # Emulate - if case ∈ ["RF-nosvd-nonsep", "RF-nosvd-sep"] + if case ∈ ["RF-nonsep"] encoder_schedule = [] # no encoding else encoder_schedule = (decorrelate_structure_mat(), "out") @@ -224,7 +204,7 @@ function main() optimize_hyperparameters!(emulator) # diagnostics - if case ∈ ["RF-svd-nonsep", "RF-nosvd-nonsep", "RF-svd-sep"] + if case ∈ ["RF-nonsep" "RF-lr-lr"] push!(opt_diagnostics, get_optimizer(mlt)[1]) #length-1 vec of vec -> vec end diff --git a/examples/Emulator/Regression_2d_2d/Project.toml b/examples/Emulator/Regression_2d_2d/Project.toml index d1ef6639c..149f5226a 100644 --- a/examples/Emulator/Regression_2d_2d/Project.toml +++ b/examples/Emulator/Regression_2d_2d/Project.toml @@ -1,8 +1,10 @@ [deps] CalibrateEmulateSample = "95e48a1f-0bec-4818-9538-3db4340308e3" +ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" +JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" diff --git a/examples/Emulator/Regression_2d_2d/compare_regression.jl b/examples/Emulator/Regression_2d_2d/compare_regression.jl index 116515d45..2becbeeea 100644 --- a/examples/Emulator/Regression_2d_2d/compare_regression.jl +++ b/examples/Emulator/Regression_2d_2d/compare_regression.jl @@ -32,8 +32,7 @@ end function main() - rng_seed = 41 - Random.seed!(rng_seed) + rng = MersenneTwister(41) output_directory = joinpath(@__DIR__, "output") if !isdir(output_directory) mkdir(output_directory) @@ -47,7 +46,8 @@ function main() "rf-nonsep", ] - encoder_names = ["no-proc", "sample-struct-proc", "sample-proc", "struct-mat-proc", "combined-proc", "cca-proc"] + encoder_names = + ["no-proc", "sample-struct-proc", "sample-proc", "struct-mat-proc", "combined-proc", "cca-proc", "minmax"] encoders = [ [], # no proc nothing, # default proc @@ -55,24 +55,35 @@ function main() (decorrelate_structure_mat(), "in_and_out"), (decorrelate(), "in_and_out"), (canonical_correlation(), "in_and_out"), + (minmax_scale(), "in_and_out"), ] # USER CHOICES - encoder_mask = [3, 4, 6] # best performing - case_mask = [1, 3, 4, 5] # (KEEP set to 1:length(cases) when pushing for buildkite) + encoder_mask = [1, 7] # best performing + case_mask = [1, 3, 5] + # scales the problem to different in-out domains (but g gives the same "function" on all domain sizes + out_scaling = 1e1 # scales range(g1) to out_scaling + in_scaling = 1.0 # scales range(x) to in_scaling + + # (necessary to get the ranges correct above) + out_scaling /= 40.0 + in_scaling /= 4.0 * pi + #problem n = 200 # number of training points p = 2 # input dim d = 2 # output dim - prior_cov = (pi^2) * I(p) - X = rand(MvNormal(zeros(p), prior_cov), n) + prior_cov = (pi^2) * I(p) * in_scaling * in_scaling + X = rand(rng, MvNormal(zeros(p), prior_cov), n) + + @info "problem scaling " in_scaling out_scaling # G(x1, x2) - g1(x) = sin.(x[1, :]) .+ cos.(2 * x[2, :]) - g2(x) = 2 * sin.(2 * x[1, :]) .- 3 * cos.(x[2, :]) - g1(x, y) = sin(x) + 2 * cos(2 * y) - g2(x, y) = 2 * sin(2 * x) - 3 * cos(y) + g1(x) = out_scaling * 10 * (sin.(1 / in_scaling * x[1, :]) .+ cos.(1 / in_scaling * 2 * x[2, :])) + g2(x) = out_scaling * 10 * (2 * sin.(1 / in_scaling * 2 * x[1, :]) .- 3 * cos.(1 / in_scaling * x[2, :])) + g1(x, y) = out_scaling * 10 * (sin(1 / in_scaling * x) + 2 * cos(1 / in_scaling * 2 * y)) + g2(x, y) = out_scaling * 10 * (2 * sin(1 / in_scaling * 2 * x) - 3 * cos(1 / in_scaling * y)) g1x = g1(X) g2x = g2(X) gx = zeros(2, n) @@ -81,8 +92,10 @@ function main() # Add noise η μ = zeros(d) - Σ = 0.1 * [[0.4, -0.3] [-0.3, 0.5]] # d x d - noise_samples = rand(MvNormal(μ, Σ), n) + #Σ = 0.1 * [[0.4, -0.3] [-0.3, 0.5]] # d x d + #Σ = 0.001 * [[0.4, -0.0] [-0.0, 0.5]] # d x d + Σ = out_scaling * out_scaling * 1.0 * I + noise_samples = rand(rng, MvNormal(μ, Σ), n) # y = G(x) + η Y = gx .+ noise_samples @@ -93,8 +106,8 @@ function main() #plot training data with and without noise if plot_flag n_pts = 200 - x1 = range(-2 * pi, stop = 2 * π, length = n_pts) - x2 = range(-2 * pi, stop = 2 * π, length = n_pts) + x1 = range(-2 * pi * in_scaling, stop = 2 * π * in_scaling, length = n_pts) + x2 = range(-2 * pi * in_scaling, stop = 2 * π * in_scaling, length = n_pts) p1 = contourf(x1, x2, g1.(x1', x2), c = :cividis, xlabel = "x1", ylabel = "x2", aspect_ratio = :equal) scatter!(p1, X[1, :], X[2, :], c = :black, ms = 3, label = "train point") @@ -123,10 +136,16 @@ function main() # common Gaussian feature setup pred_type = YType() - # common random feature setup - n_features = 200 - optimizer_options = - Dict("n_iteration" => 5, "n_features_opt" => 200, "n_ensemble" => 80, "cov_sample_multiplier" => 5.0) + # common random feature setup #large n_features, small n_features_opt + n_features = 500 + optimizer_options = Dict( + "n_iteration" => 15, + "n_features_opt" => 100, + "n_ensemble" => 100, + "cov_sample_multiplier" => 10.0, + "verbose" => true, + "cov_correction" => "nice", + ) nugget = 1e-8 # data processing schedule @@ -144,7 +163,8 @@ function main() mlt = ScalarRandomFeatureInterface( n_features, p, - kernel_structure = SeparableKernel(LowRankFactor(2, nugget), OneDimFactor()), + kernel_structure = SeparableKernel(LowRankFactor(1, nugget), OneDimFactor()), + #kernel_structure = SeparableKernel(DiagonalFactor(nugget), OneDimFactor()), optimizer_options = optimizer_options, ) elseif case == "rf-lr-lr" @@ -152,7 +172,7 @@ function main() n_features, p, d, - kernel_structure = SeparableKernel(LowRankFactor(2, nugget), LowRankFactor(2, nugget)), + kernel_structure = SeparableKernel(DiagonalFactor(nugget), LowRankFactor(2, nugget)), optimizer_options = deepcopy(optimizer_options), ) elseif case == "rf-nonsep" @@ -160,7 +180,7 @@ function main() n_features, p, d, - kernel_structure = NonseparableKernel(LowRankFactor(4, nugget)), + kernel_structure = NonseparableKernel(LowRankFactor(1, nugget)), optimizer_options = deepcopy(optimizer_options), ) @@ -179,15 +199,21 @@ function main() # Plot mean and variance of the predicted observables y1 and y2 # For this, we generate test points on a x1-x2 grid. n_pts = 200 - x1 = range(-2 * π, stop = 2 * π, length = n_pts) - x2 = range(-2 * π, stop = 2 * π, length = n_pts) + x1 = range(-2 * π * in_scaling, stop = 2 * π * in_scaling, length = n_pts) + x2 = range(-2 * π * in_scaling, stop = 2 * π * in_scaling, length = n_pts) X1, X2 = meshgrid(x1, x2) + # Input for predict has to be of size input_dim x N_samples inputs = permutedims(hcat(X1[:], X2[:]), (2, 1)) em_mean, em_cov = predict(emulator, inputs, transform_to_real = true) println("end predictions at ", n_pts * n_pts, " points") + g1_true = g1(inputs) + g1_true_grid = reshape(g1_true, n_pts, n_pts) + g2_true = g2(inputs) + g2_true_grid = reshape(g2_true, n_pts, n_pts) + g_true_grids = [g1_true_grid, g2_true_grid] #plot predictions for y_i in 1:d diff --git a/examples/GCM/emulate_sample_script.jl b/examples/GCM/emulate_sample_script.jl index 0c53cf4a5..370996d15 100644 --- a/examples/GCM/emulate_sample_script.jl +++ b/examples/GCM/emulate_sample_script.jl @@ -81,7 +81,6 @@ function main() # setup random features eki_options_override = Dict( "verbose" => true, - "tikhonov" => 0.0, "scheduler" => DataMisfitController(on_terminate = "continue"), "n_iteration" => 10, "multithread" => "ensemble", diff --git a/examples/Lorenz/calibrate_spatial_dep.jl b/examples/Lorenz/calibrate_spatial_dep.jl index 25152046c..a9801db84 100644 --- a/examples/Lorenz/calibrate_spatial_dep.jl +++ b/examples/Lorenz/calibrate_spatial_dep.jl @@ -169,12 +169,12 @@ x0 = x_spun_up[:, end] #last element of the run is the initial condition for cr #Creating sythetic data -T = 50.0 +T = 24.0 ny = nx * 2 #number of data points lorenz_config_settings = LorenzConfig(dt, T) # construct how we compute Observations -T_start = T - 25.0 +T_start = T - 20.0 T_end = T observation_config = ObservationConfig(T_start, T_end) @@ -196,9 +196,10 @@ y_ens = hcat( obs_noise_cov = cov(y_ens, dims = 2) + 1e-2 * I y = y_ens[:, 1] +#Prior covariance + pl = 4.0 psig = 5.0 -#Prior covariance B = zeros(nx, nx) for ii in 1:nx for jj in 1:nx @@ -207,6 +208,12 @@ for ii in 1:nx end B_sqrt = sqrt(B) +#= +psig = 5.0 +B = psig^2*I +B_sqrt = sqrt(B) +=# + #Prior mean mu = 4.0 * ones(nx) diff --git a/examples/Lorenz/emulate_sample.jl b/examples/Lorenz/emulate_sample.jl index b204e1df9..1f13fc177 100644 --- a/examples/Lorenz/emulate_sample.jl +++ b/examples/Lorenz/emulate_sample.jl @@ -39,6 +39,7 @@ function main() "scheduler" => DataMisfitController(terminate_at = 100.0), "cov_sample_multiplier" => 1.0, "n_iteration" => 20, + "n_features_opt" => 100, ) # we do not want termination, as our priors have relatively little interpretation diff --git a/examples/Lorenz/emulate_sample_spatial_dep.jl b/examples/Lorenz/emulate_sample_spatial_dep.jl index 04178f72a..9cffec89e 100644 --- a/examples/Lorenz/emulate_sample_spatial_dep.jl +++ b/examples/Lorenz/emulate_sample_spatial_dep.jl @@ -25,16 +25,18 @@ function main() cases = [ "GP", "RF-scalar", # diagonalize, train scalar RF, don't asume diag inputs + "RF-nonsep", # diagonalize, train scalar RF, don't asume diag inputs ] #### CHOOSE YOUR CASE: - mask = [1]# 1:1 # e.g. 1:2 or [2] + mask = [2] # e.g. 1:2 or [3] for (case) in cases[mask] println("case: ", case) min_iter = 1 - max_iter = 7 # number of EKP iterations to use data from is at most this + skip_iter = 1 + max_iter = 6 # number of EKP iterations to use data from is at most this #### @@ -71,7 +73,6 @@ function main() # Emulate-sample settings # choice of machine-learning tool in the emulation stage - nugget = 1e-3 if case == "GP" gppackage = Emulators.GPJL() pred_type = Emulators.YType() @@ -82,15 +83,17 @@ function main() noise_learn = false, ) elseif case == "RF-scalar" + nugget = 1e-6 overrides = Dict( - # "verbose" => true, + "verbose" => true, "scheduler" => DataMisfitController(terminate_at = 1000.0), "cov_sample_multiplier" => 1.0, - "n_iteration" => 8, - "n_features_opt" => 40, + "n_iteration" => 10, + "n_features_opt" => 160, ) - n_features = 80 - kernel_structure = SeparableKernel(LowRankFactor(1, nugget), OneDimFactor()) + n_features = 200 + # kernel_structure = SeparableKernel(LowRankFactor(1,nugget), OneDimFactor()) + kernel_structure = SeparableKernel(DiagonalFactor(nugget), OneDimFactor()) mlt = ScalarRandomFeatureInterface( n_features, n_params, @@ -98,6 +101,27 @@ function main() kernel_structure = kernel_structure, optimizer_options = overrides, ) + elseif case ∈ ["RF-nonsep"] + nugget = 1e-6 + overrides = Dict( + "verbose" => true, + "scheduler" => DataMisfitController(terminate_at = 1000.0), + "cov_sample_multiplier" => 1.0, + "n_iteration" => 10, + "n_features_opt" => 160, + ) + kernel_structure = NonseparableKernel(LowRankFactor(1, nugget)) + n_features = 200 + + mlt = VectorRandomFeatureInterface( + n_features, + n_params, + output_dim, + rng = rng, + kernel_structure = kernel_structure, + optimizer_options = overrides, + ) + else throw(ArgumentError("case $(case) not recognised, please choose from the list $(cases)")) end @@ -106,8 +130,10 @@ function main() # Get training points from the EKP iteration number in the second input term N_iter = min(max_iter, length(get_u(ekpobj)) - 1) # number of paired iterations taken from EKP min_iter = min(max_iter, max(1, min_iter)) - input_output_pairs = Utilities.get_training_points(ekpobj, min_iter:(N_iter - 1)) + input_output_pairs = Utilities.get_training_points(ekpobj, min_iter:skip_iter:(N_iter - 1)) input_output_pairs_test = Utilities.get_training_points(ekpobj, N_iter:(length(get_u(ekpobj)) - 1)) # "next" iterations + @info "Train iterations: $(min_iter:skip_iter:(N_iter - 1))" + @info "Test iterations: $(N_iter:(length(get_u(ekpobj)) - 1))" # Save data @save joinpath(data_save_directory, "input_output_pairs.jld2") input_output_pairs @@ -144,7 +170,7 @@ function main() emulator = Emulator( mlt, input_output_pairs; - encoder_schedule = encoder_schedule, + encoder_schedule = deepcopy(encoder_schedule), encoder_kwargs = (; prior_cov = cov(priors), obs_noise_cov = Γy), ) optimize_hyperparameters!(emulator) diff --git a/examples/Lorenz/plot_spatial_dep.jl b/examples/Lorenz/plot_spatial_dep.jl index 4265469b7..f63f57e27 100644 --- a/examples/Lorenz/plot_spatial_dep.jl +++ b/examples/Lorenz/plot_spatial_dep.jl @@ -1,8 +1,8 @@ # some context values (from calibrate_spatial_dep.jl) # some context calues from emulate_sample_spatial_dep.jl -cases = ["GP", "RF-scalar"] -case = cases[1] +cases = ["GP", "RF-scalar", "RF-nonsep"] +case = cases[2] # load packages # CES @@ -19,6 +19,7 @@ using CalibrateEmulateSample.EnsembleKalmanProcesses using CalibrateEmulateSample.EnsembleKalmanProcesses.ParameterDistributions const EKP = EnsembleKalmanProcesses + # load homedir = pwd() println(homedir) @@ -29,6 +30,8 @@ truth_params = load(data_file)["truth_params_constrained"] final_params = load(data_file)["final_params_constrained"] posterior = load(data_file)["posterior"] +prior_file = joinpath(data_save_directory, "priors_spatial_dep.jld2") +prior = load(prior_file)["prior"] ekp_file = joinpath(data_save_directory, "ekp_spatial_dep.jld2") ekpobj = load(ekp_file)["ekpobj"] @@ -49,9 +52,9 @@ quantiles = reduce(hcat, [quantile(row, [0.05, 0.5, 0.95]) for row in eachrow(co gr(size = (2 * 1.6 * 600, 600), guidefontsize = 18, tickfontsize = 16, legendfontsize = 16) p1 = plot( range(0, nx - 1, step = 1), - [truth_params final_params], - label = ["solution" "EKI"], - color = [:black :lightgreen], + truth_params, + label = "solution", + color = :black, linewidth = 4, xlabel = "Spatial index", ylabel = "Forcing (input)", @@ -72,7 +75,18 @@ p2 = plot( bottom_margin = 15mm, xticks = (Int.(0:10:ny), [0, 10, 20, 30, (40, 0), 10, 20, 30, 40]), ) -plot!(p2, 1:length(y), get_g_mean_final(ekpobj), label = "mean-final-output", color = :lightgreen, linewidth = 4) + +l = @layout [a b] +plt = plot(p1, p2, layout = l) + +savefig(plt, figure_save_directory * "data_spatial_dep.png") +savefig(plt, figure_save_directory * "data_spatial_dep.pdf") + +p3 = deepcopy(p1) +p4 = deepcopy(p2) + +plot!(p1, range(0, nx - 1, step = 1), final_params, label = "mean ensemble input", color = :lightgreen, linewidth = 4) +plot!(p2, 1:length(y), get_g_mean_final(ekpobj), label = "mean ensemble output", color = :lightgreen, linewidth = 4) l = @layout [a b] plt = plot(p1, p2, layout = l) @@ -80,6 +94,64 @@ plt = plot(p1, p2, layout = l) savefig(plt, figure_save_directory * "solution_spatial_dep_ens$(N_ens).png") savefig(plt, figure_save_directory * "solution_spatial_dep_ens$(N_ens).pdf") +# +plot!( + p3, + range(0, nx - 1, step = 1), + get_ϕ_final(prior, ekpobj)[:, 1], + label = "ensemble inputs", + color = :lightgreen, + linewidth = 4, + linealpha = 0.1, +) + +plot!( + p3, + range(0, nx - 1, step = 1), + get_ϕ_final(prior, ekpobj)[:, 2:end], + label = "", + color = :lightgreen, + linewidth = 4, + linealpha = 0.1, +) +plot!( + p3, + range(0, nx - 1, step = 1), + get_ϕ(prior, ekpobj, 1)[:, 1], + label = "", + color = :lightgreen, + linewidth = 4, + linealpha = 0.1, +) + +plot!( + p3, + range(0, nx - 1, step = 1), + get_ϕ(prior, ekpobj, 1)[:, 2:end], + label = "", + color = :lightgreen, + linewidth = 4, + linealpha = 0.1, +) + +plot!( + p4, + 1:length(y), + get_g_final(ekpobj)[:, 1], + color = :lightgreen, + label = "ensemble outputs", + linewidth = 4, + linealpha = 0.1, +) +plot!(p4, 1:length(y), get_g_final(ekpobj)[:, 2:end], color = :lightgreen, label = "", linewidth = 4, linealpha = 0.1) +plot!(p4, 1:length(y), get_g(ekpobj, 1)[:, 1], color = :lightgreen, label = "", linewidth = 4, linealpha = 0.1) +plot!(p4, 1:length(y), get_g(ekpobj, 1)[:, 2:end], color = :lightgreen, label = "", linewidth = 4, linealpha = 0.1) + + +plt = plot(p3, p4, layout = l) + +savefig(plt, figure_save_directory * "solution_spatial_dep_full_ens$(N_ens).png") +savefig(plt, figure_save_directory * "solution_spatial_dep_full_ens$(N_ens).pdf") # plot - UQ results @@ -107,7 +179,86 @@ plot!( fillalpha = 0.1, ) - figpath = joinpath(figure_save_directory, "posterior_ribbons_" * case) savefig(figpath * ".png") savefig(figpath * ".pdf") + +########################## +#cases = ["GP", "RF-scalar", "RF-nonsep"] + +# load +homedir = pwd() +println(homedir) +data_file_GP = joinpath(data_save_directory, "posterior_$(cases[1]).jld2") +data_file_RF = joinpath(data_save_directory, "posterior_$(case).jld2") +truth_params_GP = load(data_file_GP)["truth_params_constrained"] +final_params_GP = load(data_file_GP)["final_params_constrained"] +posterior_GP = load(data_file_GP)["posterior"] +truth_params_RF = load(data_file_RF)["truth_params_constrained"] +final_params_RF = load(data_file_RF)["final_params_constrained"] +posterior_RF = load(data_file_RF)["posterior"] + +ekp_file = joinpath(data_save_directory, "ekp_spatial_dep.jld2") +ekpobj = load(ekp_file)["ekpobj"] +N_ens = get_N_ens(ekpobj) +nx = length(truth_params) +y = get_obs(ekpobj) +ny = length(y) +# get samples and quantiles from posterior +param_names_GP = get_name(posterior_GP) +posterior_samples_GP = vcat([get_distribution(posterior_GP)[name] for name in get_name(posterior_GP)]...) #samples are columns +constrained_posterior_samples_GP = + mapslices(x -> transform_unconstrained_to_constrained(posterior_GP, x), posterior_samples_GP, dims = 1) + +quantiles_GP = reduce(hcat, [quantile(row, [0.05, 0.5, 0.95]) for row in eachrow(constrained_posterior_samples_GP)])' # rows are quantiles for row of posterior samples + +param_names_RF = get_name(posterior_RF) +posterior_samples_RF = vcat([get_distribution(posterior_RF)[name] for name in get_name(posterior_RF)]...) #samples are columns +constrained_posterior_samples_RF = + mapslices(x -> transform_unconstrained_to_constrained(posterior_RF, x), posterior_samples_RF, dims = 1) + +quantiles_RF = reduce(hcat, [quantile(row, [0.05, 0.5, 0.95]) for row in eachrow(constrained_posterior_samples_RF)])' # rows are quantiles for row of posterior samples + +# plot - UQ results - both + +gr(size = (1.6 * 600, 600), guidefontsize = 18, tickfontsize = 16, legendfontsize = 16) +p1 = plot( + range(0, nx - 1, step = 1), + [truth_params final_params], + label = ["solution" "EKI-opt"], + color = [:black :lightgreen], + linewidth = 4, + xlabel = "Spatial index", + ylabel = "Forcing (input)", + left_margin = 15mm, + bottom_margin = 15mm, +) + +plot!( + p1, + range(0, nx - 1, step = 1), + quantiles_GP[:, 2], # median of all vals + color = :blue, + label = "GP posterior", + linewidth = 4, + ribbon = [quantiles_GP[:, 2] - quantiles_GP[:, 1] quantiles_GP[:, 3] - quantiles_GP[:, 2]], + linealpha = 0.5, + fillalpha = 0.1, +) + +plot!( + p1, + range(0, nx - 1, step = 1), + quantiles_RF[:, 2], # median of all vals + color = :red, + label = "posterior_RF", + linewidth = 4, + ribbon = [quantiles_RF[:, 2] - quantiles_RF[:, 1] quantiles_RF[:, 3] - quantiles_RF[:, 2]], + linealpha = 0.5, + fillalpha = 0.1, +) + + +figpath = joinpath(figure_save_directory, "posterior_ribbons_$(cases[1])_$(case)") +savefig(figpath * ".png") +savefig(figpath * ".pdf") diff --git a/src/RandomFeature.jl b/src/RandomFeature.jl index 7a07eae8a..4383d7f56 100644 --- a/src/RandomFeature.jl +++ b/src/RandomFeature.jl @@ -165,7 +165,7 @@ calculates the number of hyperparameters generated by the choice of covariance s calculate_n_hyperparameters(d::Int, odf::OneDimFactor) = 0 calculate_n_hyperparameters(d::Int, df::DiagonalFactor) = d calculate_n_hyperparameters(d::Int, cf::CholeskyFactor) = Int(d * (d + 1) / 2) -calculate_n_hyperparameters(d::Int, lrf::LowRankFactor) = Int(rank(lrf) + d * rank(lrf)) +calculate_n_hyperparameters(d::Int, lrf::LowRankFactor) = Int(d * rank(lrf) + d) calculate_n_hyperparameters(d::Int, hlrf::HierarchicalLowRankFactor) = Int(rank(hlrf) * (rank(hlrf) + 1) / 2 + d * rank(hlrf)) @@ -208,12 +208,11 @@ end function hyperparameters_from_flat(x::V, lrf::LowRankFactor) where {V <: AbstractVector} r = rank(lrf) - d = Int(length(x) / r - 1) # l = r + d*r + d = Int(length(x) / (r + 1)) # l = d + d*r = (r+1)*d - D = Diagonal(x[1:r]) # D acts like sing.vals.^2 - U = reshape(x[(r + 1):end], d, r) - qrU = qr(U) # qrU.Q is semi-orthogonal in first r columns, i.e. Q^T.Q = I - return qrU.Q[:, 1:r] * D * qrU.Q[:, 1:r]' + get_eps(lrf) * I + D = Diagonal(x[1:d]) # D acts like sing.vals.^2 + U = reshape(x[(d + 1):end], d, r) + return U * U' + D + get_eps(lrf) * I end @@ -236,46 +235,52 @@ function hyperparameters_from_flat(x::V, hlrf::HierarchicalLowRankFactor) where end -""" -$(DocStringExtensions.TYPEDSIGNATURES) -builds a prior distribution for the kernel hyperparameters to initialize optimization. -""" -build_default_prior(name::SS, n_hp::Int, odf::OneDimFactor) where {SS <: AbstractString} = nothing +build_default_prior(name::SS, n_hp::Int, odf::OneDimFactor; kwargs...) where {SS <: AbstractString} = nothing function build_default_prior(name::SS, n_hp::Int, df::DiagonalFactor) where {SS <: AbstractString} return ParameterDistribution( Dict( "name" => "$(name)_diagonal", - "distribution" => VectorOfParameterized(repeat([Normal(0, 3)], n_hp)), - "constraint" => repeat([bounded_below(0.0)], n_hp), + "distribution" => VectorOfParameterized(repeat([Normal(log(sqrt(0.1 / n_hp)), 2)], n_hp)), + "constraint" => repeat([bounded_below(1e-8 / n_hp)], n_hp), ), ) end function build_default_prior(name::SS, n_hp::Int, df::CholeskyFactor) where {SS <: AbstractString} - return constrained_gaussian("$(name)_cholesky", 0.0, 5.0, -Inf, Inf, repeats = n_hp) + return constrained_gaussian("$(name)_cholesky", 0.0, 0.1 / n_hp, -Inf, Inf, repeats = n_hp) end function build_default_prior(name::SS, n_hp::Int, lrf::LowRankFactor) where {SS <: AbstractString} r = rank(lrf) - d = Int(n_hp / r - 1) + d = Int(n_hp / (r + 1)) # n_hp = d + d*r D = ParameterDistribution( Dict( "name" => "$(name)_lowrank_diagonal", - "distribution" => VectorOfParameterized(repeat([Normal(0, 3)], r)), - "constraint" => repeat([bounded_below(0.0)], r), + "distribution" => VectorOfParameterized(repeat([Normal(log(sqrt(0.1 / d)), 2)], d)), + # "distribution" => VectorOfParameterized(repeat([Normal(log(1.0 / d), 2)], d)), + "constraint" => repeat([bounded_below(1e-8 / d)], d), ), ) - U = constrained_gaussian("$(name)_lowrank_U", 0.0, 100.0, -Inf, Inf, repeats = Int(d * r)) + # U = constrained_gaussian("$(name)_lowrank_U", 0.0, 10.0 / (d * r), -Inf, Inf, repeats = Int(d * r)) + U = constrained_gaussian("$(name)_lowrank_U", 0.0, sqrt(0.1 / d) / r, -Inf, Inf, repeats = Int(d * r)) + # D scales with sqrt(1/d), U scales with 1/r,sqrt(1/d) return combine_distributions([D, U]) end function build_default_prior(name::SS, n_hp::Int, hlrf::HierarchicalLowRankFactor) where {SS <: AbstractString} r = rank(hlrf) d = Int(n_hp / r - (r + 1) / 2) - L = constrained_gaussian("$(name)_lowrank_Kchol", 0.0, 5.0, -Inf, Inf, repeats = Int(r * (r + 1) / 2)) - U = constrained_gaussian("$(name)_lowrank_U", 0.0, 100.0, -Inf, Inf, repeats = Int(d * r)) + L = constrained_gaussian( + "$(name)_lowrank_Kchol", + 0.0, + 1.0 / (r * (r + 1) / 2), + -Inf, + Inf, + repeats = Int(r * (r + 1) / 2), + ) + U = constrained_gaussian("$(name)_lowrank_U", 0.0, sqrt(0.1 / d) / r, -Inf, Inf, repeats = Int(d * r)) return combine_distributions([L, U]) end @@ -290,14 +295,13 @@ function calculate_n_hyperparameters( n_hp_in = calculate_n_hyperparameters(input_dim, input_cov_structure) output_cov_structure = get_output_cov_structure(kernel_structure) n_hp_out = calculate_n_hyperparameters(output_dim, output_cov_structure) - n_scalings = 1 # if both in and output dim are "OneDimFactors" we still learn 1 hp if n_hp_in + n_hp_out == 0 - return 1 + n_scalings + return 1 end - return n_hp_in + n_hp_out + n_scalings + return n_hp_in + n_hp_out end function calculate_n_hyperparameters(input_dim::Int, kernel_structure::KST) where {KST <: KernelStructureType} @@ -312,18 +316,22 @@ function calculate_n_hyperparameters( ) where {NK <: NonseparableKernel} cov_structure = get_cov_structure(kernel_structure) n_hp = calculate_n_hyperparameters(input_dim * output_dim, cov_structure) - n_scalings = 1 + # if both in and output dim are "OneDimFactors" we still learn 1 hp if n_hp == 0 - return 1 + n_scalings + return 1 end - return n_hp + n_scalings + return n_hp end +""" +$(DocStringExtensions.TYPEDSIGNATURES) +Builds a prior distribution for the kernel hyperparameters to initialize optimization. The parameter distributions built from these priors will be scaled such that the input and output range of the data is O(1). +""" function build_default_prior(input_dim::Int, output_dim::Int, kernel_structure::SK) where {SK <: SeparableKernel} input_cov_structure = get_input_cov_structure(kernel_structure) output_cov_structure = get_output_cov_structure(kernel_structure) @@ -331,21 +339,20 @@ function build_default_prior(input_dim::Int, output_dim::Int, kernel_structure:: n_hp_in = calculate_n_hyperparameters(input_dim, input_cov_structure) input_prior = build_default_prior("input", n_hp_in, input_cov_structure) n_hp_out = calculate_n_hyperparameters(output_dim, output_cov_structure) - output_prior = build_default_prior("output", n_hp_out, output_cov_structure) - - scaling_kern = constrained_gaussian("sigma", 1, 5, 0, Inf) + output_prior = build_default_prior("output", n_hp_out, output_cov_structure) # use output scale only for sigma # We only use OneDimFactor values, default to input if both in and output dimension are one dimensional.Otherwise they are ignored if isnothing(input_prior) && isnothing(output_prior) onedim_prior = constrained_gaussian("onedim", 1.0, 1.0, 0.0, Inf) - return combine_distributions([onedim_prior, scaling_kern]) + return onedim_prior elseif isnothing(input_prior) - return combine_distributions([output_prior, scaling_kern]) + return output_prior elseif isnothing(output_prior) - return combine_distributions([input_prior, scaling_kern]) + return input_prior else - return combine_distributions([input_prior, output_prior, scaling_kern]) + return combine_distributions([input_prior, output_prior]) end + end function build_default_prior(input_dim::Int, kernel_structure::KST) where {KST <: KernelStructureType} @@ -362,8 +369,8 @@ function build_default_prior(input_dim::Int, output_dim::Int, kernel_structure:: pd_kern = build_default_prior("full", n_hp, cov_structure) end - scaling_kern = constrained_gaussian("sigma", 1, 5, 0, Inf) - return combine_distributions([pd_kern, scaling_kern]) + return pd_kern + end function hyperparameters_from_flat( @@ -440,6 +447,8 @@ function calculate_mean_cov_and_coeffs( cov_store::A, buffer::A, multithread_type::MT, + prior_in_scale, + prior_out_scale, ) where { RFI <: RandomFeatureInterface, RNG <: AbstractRNG, @@ -473,6 +482,8 @@ function calculate_mean_cov_and_coeffs( input_dim, output_dim, multithread_type, + prior_in_scale, + prior_out_scale, ) fitted_features = RF.Methods.fit(rfm, io_train_cost, decomposition_type = decomp_type) @@ -487,15 +498,11 @@ function calculate_mean_cov_and_coeffs( buffer, tullio_threading = thread_opt, ) - # sizes (output_dim x n_test), (output_dim x output_dim x n_test) ## TODO - the theory states that the following should be set: scaled_coeffs = sqrt(1 / (n_features)) * RF.Methods.get_coeffs(fitted_features) - #scaled_coeffs = 1e-3 * rand(n_features)#overwrite with noise... - # However the convergence is much improved with setting this to zero: - #scaled_coeffs = 0 if decomp_type == "cholesky" chol_fac = RF.Methods.get_decomposition(RF.Methods.get_feature_factors(fitted_features)).L @@ -505,18 +512,17 @@ function calculate_mean_cov_and_coeffs( complexity = sum(log, svd_singval) # note this is log(abs(det)) end complexity = sqrt(complexity) + return scaled_coeffs, complexity end - - """ $(DocStringExtensions.TYPEDSIGNATURES) Calculate the empirical covariance, additionally applying a shrinkage operator (here the Ledoit Wolf 2004 shrinkage operation). Known to have better stability properties than Monte-Carlo for low sample sizes """ -function shrinkage_cov(sample_mat::AA) where {AA <: AbstractMatrix} +function shrinkage_cov(sample_mat::AA; cov_or_corr = "cov", verbose = false) where {AA <: AbstractMatrix} n_out, n_sample = size(sample_mat) # de-mean (as we will use the samples directly for calculation of β) @@ -524,22 +530,42 @@ function shrinkage_cov(sample_mat::AA) where {AA <: AbstractMatrix} # Ledoit Wolf shrinkage to I # get sample covariance - Γ = cov(sample_mat_zeromean, dims = 2) + if cov_or_corr == "cov" + Γ = cov(sample_mat_zeromean, dims = 2) + ss = sample_mat_zeromean + else # "corr" + Γcov = cov(sample_mat_zeromean, dims = 2) + bd_tol = 1e8 * eps() + v = sqrt.(diag(Γcov)) + V = Diagonal(v) #stds + V_inv = inv(V) + Γ = clamp.(V_inv * Γcov * V_inv, -1 + bd_tol, 1 - bd_tol) # full corr + ss = V_inv * sample_mat_zeromean + + end # estimate opt shrinkage μ_shrink = 1 / n_out * tr(Γ) δ_shrink = norm(Γ - μ_shrink * I)^2 / n_out # (scaled) frob norm of Γ_m #once de-meaning, we need to correct the sample covariance with an n_sample -> n_sample-1 - β_shrink = sum([norm(c * c' - Γ)^2 / n_out for c in eachcol(sample_mat_zeromean)]) / (n_sample - 1)^2 + β_shrink = sum([norm(c * c' - Γ)^2 / n_out for c in eachcol(ss)]) / (n_sample - 1)^2 γ_shrink = min(β_shrink / δ_shrink, 1) # clipping is typically rare + if β_shrink / δ_shrink > 1 && verbose + @warn "clipping applied; forcing shrinkage parameters β/δ=$(β_shrink/δ_shrink) -> 1." + end # γμI + (1-γ)Γ Γ .*= (1 - γ_shrink) for i in 1:n_out Γ[i, i] += γ_shrink * μ_shrink end - - @info "Shrinkage scale: $(γ_shrink), (0 = none, 1 = revert to scaled Identity)\n shrinkage covariance condition number: $(cond(Γ))" - return Γ + if verbose + @info "Shrinkage scale: $(γ_shrink), (0 = none, 1 = revert to scaled Identity)\n shrinkage covariance condition number: $(cond(Γ))" + end + if cov_or_corr == "cov" + return Γ + else + return V * Γ * V + end end @@ -548,7 +574,7 @@ $(DocStringExtensions.TYPEDSIGNATURES) Calculate the empirical covariance, additionally applying the Noise Informed Covariance Estimator (NICE) Vishnay et al. 2024. """ -function nice_cov(sample_mat::AA, δ::FT = 1.0) where {AA <: AbstractMatrix, FT <: Real} +function nice_cov(sample_mat::AA; δ::FT = 1.0, verbose = false) where {AA <: AbstractMatrix, FT <: Real} n_sample_cov = size(sample_mat, 2) Γ = cov(sample_mat, dims = 2) @@ -587,12 +613,37 @@ function nice_cov(sample_mat::AA, δ::FT = 1.0) where {AA <: AbstractMatrix, FT end end out = posdef_correct(V * corr * V) # rebuild the cov matrix - @info "NICE-adjusted covariance condition number: $(cond(out))" + if verbose + @info "NICE-adjusted covariance condition number: $(cond(out))" + end return out end - +""" +Build covariance from the three blocks, blockmeans, coeffl2norm and complexity with additional regularization +""" +function correct_covariance( + blockmeans::AM, + coeffl2norm::AM, + complexity::AM, + cov_correction::AS; + verbose = false, +) where {AM <: AbstractMatrix, AS <: AbstractString} + if cov_correction == "shrinkage" + sample_mat = vcat(blockmeans, coeffl2norm, complexity) + return shrinkage_cov(sample_mat; verbose = verbose) + elseif cov_correction == "shrinkage_corr" + sample_mat = vcat(blockmeans, coeffl2norm, complexity) + return shrinkage_cov(sample_mat; cov_or_corr = "corr", verbose = verbose) + elseif cov_correction == "nice" + sample_mat = vcat(blockmeans, coeffl2norm, complexity) + return nice_cov(sample_mat; verbose = verbose) + else + sample_mat = vcat(blockmeans, coeffl2norm, complexity) + return cov(sample_mat, dims = 2) + end +end """ $(DocStringExtensions.TYPEDSIGNATURES) @@ -611,9 +662,12 @@ function estimate_mean_and_coeffnorm_covariance( io_pairs::PairedDataContainer, n_samples::Int, decomp_type::S, - multithread_type::TullioThreading; + multithread_type::TullioThreading, + prior_in_scale, + prior_out_scale; repeats::Int = 1, cov_correction = "shrinkage", + verbose = false, ) where { RFI <: RandomFeatureInterface, RNG <: AbstractRNG, @@ -633,7 +687,9 @@ function estimate_mean_and_coeffnorm_covariance( buffer = zeros(n_test, output_dim, n_features) complexity = zeros(1, n_samples) coeffl2norm = zeros(1, n_samples) - println("estimate cov with " * string(n_samples) * " iterations...") + if verbose + println("estimate cov with " * string(n_samples) * " iterations...") + end for i in ProgressBar(1:n_samples) for j in 1:repeats @@ -652,6 +708,8 @@ function estimate_mean_and_coeffnorm_covariance( moc_tmp, buffer, multithread_type, + prior_in_scale, + prior_out_scale, ) # m output_dim x n_test @@ -680,15 +738,7 @@ function estimate_mean_and_coeffnorm_covariance( blockmeans[id, :] = permutedims(means[i, :, :], (2, 1)) end - sample_mat = vcat(blockmeans, coeffl2norm, complexity) - - if cov_correction == "shrinkage" - Γ = shrinkage_cov(sample_mat) - elseif cov_correction == "nice" - Γ = nice_cov(sample_mat) - else - Γ = cov(sample_mat, dims = 2) - end + Γ = correct_covariance(blockmeans, coeffl2norm, complexity, cov_correction; verbose = verbose) if !isposdef(approx_σ2) println("approx_σ2 not posdef") @@ -699,6 +749,7 @@ function estimate_mean_and_coeffnorm_covariance( end + """ $(DocStringExtensions.TYPEDSIGNATURES) @@ -715,8 +766,11 @@ function calculate_ensemble_mean_and_coeffnorm( batch_sizes::Union{Dict{S, Int}, Nothing}, io_pairs::PairedDataContainer, decomp_type::S, - multithread_type::TullioThreading; + multithread_type::TullioThreading, + prior_in_scale, + prior_out_scale; repeats::Int = 1, + verbose = false, ) where { RFI <: RandomFeatureInterface, RNG <: AbstractRNG, @@ -763,6 +817,8 @@ function calculate_ensemble_mean_and_coeffnorm( moc_tmp, buffer, multithread_type, + prior_in_scale, + prior_out_scale, ) # m output_dim x n_test # v output_dim x output_dim x n_test @@ -805,9 +861,12 @@ function estimate_mean_and_coeffnorm_covariance( io_pairs::PairedDataContainer, n_samples::Int, decomp_type::S, - multithread_type::EnsembleThreading; + multithread_type::EnsembleThreading, + prior_in_scale, + prior_out_scale; repeats::Int = 1, cov_correction = "shrinkage", + verbose = false, ) where { RFI <: RandomFeatureInterface, RNG <: AbstractRNG, @@ -819,8 +878,9 @@ function estimate_mean_and_coeffnorm_covariance( output_dim = size(get_outputs(io_pairs), 1) n_test = length(test_idx) - println("estimate cov with " * string(n_samples) * " iterations...") - + if verbose + println("estimate cov with " * string(n_samples) * " iterations...") + end nthreads = Threads.nthreads() coeffl2norm = zeros(1, n_samples) @@ -861,6 +921,8 @@ function estimate_mean_and_coeffnorm_covariance( moc_tmp[tid], buffer[tid], multithread_type, + prior_in_scale, + prior_out_scale, ) @@ -869,10 +931,9 @@ function estimate_mean_and_coeffnorm_covariance( # c n_features # cplxty 1 - # update vbles needed for cov # update vbles needed for cov means[:, i, :] += mtmp[tid] ./ repeats - coeffl2norm[1, i] += sqrt(sum(abs2, c)) / repeats + coeffl2norm[1, i] += sqrt(sum(c .^ 2)) / repeats complexity[1, i] += cplxty / repeats # update vbles needed for mean @@ -893,17 +954,7 @@ function estimate_mean_and_coeffnorm_covariance( blockmeans[id, :] = permutedims(means[i, :, :], (2, 1)) end - - sample_mat = vcat(blockmeans, coeffl2norm, complexity) - - if cov_correction == "shrinkage" - Γ = shrinkage_cov(sample_mat) - elseif cov_correction == "nice" - Γ = nice_cov(sample_mat) - else - Γ = cov(sample_mat, dims = 2) - end - + Γ = correct_covariance(blockmeans, coeffl2norm, complexity, cov_correction; verbose = verbose) if !isposdef(approx_σ2) println("approx_σ2 not posdef") approx_σ2 = posdef_correct(approx_σ2) @@ -927,8 +978,11 @@ function calculate_ensemble_mean_and_coeffnorm( batch_sizes::Union{Dict{S, Int}, Nothing}, io_pairs::PairedDataContainer, decomp_type::S, - multithread_type::EnsembleThreading; + multithread_type::EnsembleThreading, + prior_in_scale, + prior_out_scale; repeats::Int = 1, + verbose = false, ) where { RFI <: RandomFeatureInterface, RNG <: AbstractRNG, @@ -989,6 +1043,8 @@ function calculate_ensemble_mean_and_coeffnorm( moc_tmp[tid], buffer[tid], multithread_type, + prior_in_scale, + prior_out_scale, ) # m output_dim x n_test @@ -1017,7 +1073,6 @@ function calculate_ensemble_mean_and_coeffnorm( println("blockcovmat not posdef") blockcovmat = posdef_correct(blockcovmat) end - return vcat(blockmeans, coeffl2norm, complexity), blockcovmat end diff --git a/src/ScalarRandomFeature.jl b/src/ScalarRandomFeature.jl index fa59c6741..aca965e65 100644 --- a/src/ScalarRandomFeature.jl +++ b/src/ScalarRandomFeature.jl @@ -127,20 +127,19 @@ Constructs a `ScalarRandomFeatureInterface <: MachineLearningTool` interface for - `feature_decomposition = "cholesky"` - choice of how to store decompositions of random features, `cholesky` or `svd` available - `optimizer_options = nothing` - Dict of options to pass into EKI optimization of hyperparameters (defaults created in `ScalarRandomFeatureInterface` constructor): - "prior": the prior for the hyperparameter optimization - - "prior_in_scale": use this to tune the input prior scale - "n_ensemble": number of ensemble members - "n_iteration": number of eki iterations - "cov_sample_multiplier": increase for more samples to estimate covariance matrix in optimization (default 10.0, minimum 0.0) - "scheduler": Learning rate Scheduler (a.k.a. EKP timestepper) Default: DataMisfitController - - "tikhonov": tikhonov regularization parameter if >0 - "inflation": additive inflation ∈ [0,1] with 0 being no inflation - "train_fraction": e.g. 0.8 (default) means 80:20 train - test split - "n_features_opt": fix the number of features for optimization (default `n_features`, as used for prediction) - "multithread": how to multithread. "ensemble" (default) threads across ensemble members "tullio" threads random feature matrix algebra - "accelerator": use EKP accelerators (default is no acceleration) - - "verbose" => false, verbose optimizer statements - - "cov_correction" => "shrinkage", type of conditioning to improve estimated covariance (Ledoit Wolfe 03), also "nice" for (Vishny, Morzfeld et al. 2024) - - "n_cross_val_sets" => 2, train fraction creates (default 5) train-test data subsets, then use 'n_cross_val_sets' of these stacked in the loss function. If set to 0, train=test on the full data provided ignoring "train_fraction". + - "verbose" => false: verbose optimizer statements + - "cov_correction" => "nice": type of conditioning to improve estimated covariance. "shrinkage", "shrinkage_corr" (Ledoit Wolfe 03), "nice" for (Vishny, Morzfeld et al. 2024) + - "overfit" => 1.0: if > 1.0 forcibly overfit/under-regularize the optimizer cost, (vice versa for < 1.0). + - "n_cross_val_sets" => 2: train fraction creates (default 5) train-test data subsets, then use 'n_cross_val_sets' of these stacked in the loss function. If set to 0, train=test on the full data provided ignoring "train_fraction". """ function ScalarRandomFeatureInterface( n_features::Int, @@ -177,8 +176,9 @@ function ScalarRandomFeatureInterface( "verbose" => false, # verbose optimizer statements "accelerator" => EKP.NesterovAccelerator(), # acceleration with momentum "localization" => EKP.Localizers.NoLocalization(), # localization / sample error correction for small ensembles - "cov_correction" => "shrinkage", # type of conditioning to improve estimated covariance + "cov_correction" => "nice", # type of conditioning to improve estimated covariance "n_cross_val_sets" => 2, # if >1 do cross validation, else if 0 do no data splitting and no training fraction + "overfit" => 1.0, # if >1 this forcibly overfits to the data ) if !isnothing(optimizer_options) @@ -193,7 +193,9 @@ function ScalarRandomFeatureInterface( opt_tmp[key] = optimizer_opts[key] end end - @info("hyperparameter optimization with EKI configured with $opt_tmp") + if optimizer_opts["verbose"] + @info("hyperparameter optimization with EKI configured with $opt_tmp") + end return ScalarRandomFeatureInterface{S, RNG, KSType}( rfms, @@ -214,17 +216,19 @@ function hyperparameter_distribution_from_flat( x::VV, input_dim::Int, kernel_structure::SK, + prior_in_scale, ) where {VV <: AbstractVector, SK <: SeparableKernel} - M = zeros(input_dim) #scalar output U = hyperparameters_from_flat(x, input_dim, kernel_structure) - - if !isposdef(U) + # make symmetric + UU = Diagonal(vec(prior_in_scale)) * U * Diagonal(vec(prior_in_scale)) + UU = 0.5 * (UU + UU') + if !isposdef(UU) println("U not posdef - correcting") - U = posdef_correct(U) + UU = posdef_correct(UU) end - dist = MvNormal(M, U) + dist = MvNormal(zeros(input_dim), UU) pd = ParameterDistribution( Dict( "distribution" => Parameterized(dist), @@ -232,7 +236,6 @@ function hyperparameter_distribution_from_flat( "name" => "xi", ), ) - return pd end @@ -240,6 +243,7 @@ function hyperparameter_distribution_from_flat( x::VV, input_dim::Int, kernel_structure::NK, + prior_in_scale, ) where {VV <: AbstractVector, NK <: NonseparableKernel} throw( ArgumentError( @@ -261,7 +265,9 @@ function RFM_from_hyperparameters( n_features::Int, batch_sizes::Union{Dict{S, Int}, Nothing}, input_dim::Int, - multithread_type::MT; + multithread_type::MT, + prior_in_scale, + prior_out_scale; ) where { RNG <: AbstractRNG, ForVM <: Union{Real, AbstractVecOrMat}, @@ -270,14 +276,13 @@ function RFM_from_hyperparameters( MT <: MultithreadType, } - xi_hp = l[1:(end - 1)] - sigma_hp = l[end] + xi_hp = isa(l, AbstractVecOrMat) ? l[:] : [l] kernel_structure = get_kernel_structure(srfi) - pd = hyperparameter_distribution_from_flat(xi_hp, input_dim, kernel_structure) + pd = hyperparameter_distribution_from_flat(xi_hp, input_dim, kernel_structure, prior_in_scale) feature_sampler = RF.Samplers.FeatureSampler(pd, rng = rng) # Learn hyperparameters for different feature types - feature_parameters = Dict("sigma" => sigma_hp) + feature_parameters = Dict("sigma" => prior_out_scale * sqrt(2)) sff = RF.Features.ScalarFourierFeature(n_features, feature_sampler, feature_parameters = feature_parameters) thread_opt = isa(multithread_type, TullioThreading) # if we want to multithread with tullio if isnothing(batch_sizes) @@ -302,14 +307,27 @@ RFM_from_hyperparameters( batch_sizes::Union{Dict{S, Int}, Nothing}, input_dim::Int, output_dim::Int, - multithread_type::MT; + multithread_type::MT, + prior_in_scale, + prior_out_scale, ) where { RNG <: AbstractRNG, ForVM <: Union{Real, AbstractVecOrMat}, MorUS <: Union{AbstractMatrix, UniformScaling}, S <: AbstractString, MT <: MultithreadType, -} = RFM_from_hyperparameters(srfi, rng, l, regularization, n_features, batch_sizes, input_dim, multithread_type) +} = RFM_from_hyperparameters( + srfi, + rng, + l, + regularization, + n_features, + batch_sizes, + input_dim, + multithread_type, + prior_in_scale, + prior_out_scale, +) """ $(DocStringExtensions.TYPEDSIGNATURES) @@ -333,7 +351,6 @@ function build_models!( kernel_structure = get_kernel_structure(srfi) n_hp = calculate_n_hyperparameters(input_dim, kernel_structure) - rfms = get_rfms(srfi) if length(rfms) > 0 @warn "ScalarRandomFeatureInterface already built. skipping..." @@ -345,6 +362,7 @@ function build_models!( rng = get_rng(srfi) decomp_type = get_feature_decomposition(srfi) optimizer_options = get_optimizer_options(srfi) + opt_verbose_flag = optimizer_options["verbose"] optimizer = get_optimizer(srfi) # empty vector # Optimize features with EKP for each output dim @@ -376,7 +394,6 @@ function build_models!( ) end - for i in 1:n_cross_val_sets tmp = idx_shuffle[((i - 1) * n_test + 1):(i * n_test)] push!(test_idx, tmp) @@ -401,8 +418,10 @@ function build_models!( n_iteration = optimizer_options["n_iteration"] diagnostics = zeros(n_iteration, n_rfms) for i in 1:n_rfms + if opt_verbose_flag + @info "training model $i / $n_rfms" + end regularization_i = regularization[i, i] * I - io_pairs_opt = PairedDataContainer(input_values, reshape(output_values[i, :], 1, size(output_values, 2))) multithread = optimizer_options["multithread"] @@ -417,27 +436,36 @@ function build_models!( ), ) end - # [2.] Estimate covariance at mean value - prior = optimizer_options["prior"] + + # scale up the prior so that default priors are always "reasonable" + prior_in_scale = 1.0 ./ std(input_values, dims = 2) + prior_out_scale = std(output_values[i, :]) + + + prior = build_default_prior(input_dim, kernel_structure) # where prior space has changed we need to rebuild the priors if ndims(prior) > n_hp # comes from having a truncated output_dimension # TODO not really a truncation here, resetting to default - @info "Original input space of dimension $(get_input_dim(srfi)) has been truncated to $(input_dim). \n Rebuilding prior... number of hyperparameters reduced from $(ndims(prior)) to $(n_hp)." + @info "Original input space of dimension $(get_input_dim(srfi)) has been truncated to $(input_dim). \n Rebuilding RF prior... number of hyperparameters reduced from $(ndims(prior)) to $(n_hp)." prior = build_default_prior(input_dim, kernel_structure) end - + # [2.] Estimate covariance at mean value μ_hp = transform_unconstrained_to_constrained(prior, mean(prior)) cov_sample_multiplier = optimizer_options["cov_sample_multiplier"] cov_correction = optimizer_options["cov_correction"] + overfit = max(optimizer_options["overfit"], 1e-4) n_cov_samples_min = n_test + 2 n_cov_samples = Int(floor(n_cov_samples_min * max(cov_sample_multiplier, 0.0))) + println("estimating covariances with " * string(n_cov_samples) * " iterations...") + observation_vec = [] for cv_idx in 1:n_cross_val_sets + internal_Γ, approx_σ2 = estimate_mean_and_coeffnorm_covariance( srfi, rng, @@ -451,11 +479,19 @@ function build_models!( n_cov_samples, decomp_type, multithread_type, + prior_in_scale, + prior_out_scale, cov_correction = cov_correction, + verbose = opt_verbose_flag, ) - Γ = internal_Γ - Γ[1:n_test, 1:n_test] += regularization_i # + approx_σ2 + + # blocks: + Γ = deepcopy(internal_Γ) + Γ[1:n_test, 1:n_test] += regularization_i(n_test) # approx_σ2 + Γ[1:n_test, 1:n_test] /= overfit^2 # shrink the data noise artificially Γ[(n_test + 1):end, (n_test + 1):end] += I + # small features this has a larger effect - though doesn't -> I as n-> infty + if !isposdef(Γ) Γ = posdef_correct(Γ) end @@ -470,16 +506,19 @@ function build_models!( # [3.] set up EKP optimization n_ensemble = optimizer_options["n_ensemble"] n_iteration = optimizer_options["n_iteration"] - opt_verbose_flag = optimizer_options["verbose"] scheduler = optimizer_options["scheduler"] accelerator = optimizer_options["accelerator"] localization = optimizer_options["localization"] initial_params = construct_initial_ensemble(rng, prior, n_ensemble) + # bug with scalar mean o/w + prior_mean = isa(mean(prior), AbstractVector) ? mean(prior) : [mean(prior)] + prior_cov = cov(prior) + ekiobj = EKP.EnsembleKalmanProcess( initial_params, observation, - Inversion(), + TransformInversion(prior_mean, prior_cov), scheduler = scheduler, rng = rng, accelerator = accelerator, @@ -492,9 +531,9 @@ function build_models!( for i in 1:n_iteration #get parameters: - lvec = transform_unconstrained_to_constrained(prior, get_u_final(ekiobj)) + lvec = get_ϕ_final(prior, ekiobj) g_ens = zeros(n_cross_val_sets * (n_test + 2), n_ensemble) - for cv_idx in 1:n_cross_val_sets + for (iii, cv_idx) in enumerate(1:n_cross_val_sets) g_ens_tmp, _ = calculate_ensemble_mean_and_coeffnorm( srfi, @@ -508,10 +547,15 @@ function build_models!( io_pairs_opt, decomp_type, multithread_type, + prior_in_scale, + prior_out_scale, + verbose = opt_verbose_flag, ) + # useful diagnostic: + # mm = mean(g_ens_tmp, dims=2) ./ sqrt.(diag(get_obs_noise_cov(observation_vec[iii]))) + # @info mean(mm[1:end-2]) mm[end-1] mm[end] g_ens[((cv_idx - 1) * (n_test + 2) + 1):(cv_idx * (n_test + 2)), :] = g_ens_tmp end - inflation = optimizer_options["inflation"] if inflation > 0 terminated = EKP.update_ensemble!(ekiobj, g_ens, additive_inflation = true, s = inflation) # small regularizing inflation @@ -561,6 +605,8 @@ function build_models!( batch_sizes, input_dim, multithread_type, + prior_in_scale, + prior_out_scale, ) fitted_features_i = RF.Methods.fit(rfm_i, io_pairs_i, decomposition_type = decomp_type) #fit features diff --git a/src/VectorRandomFeature.jl b/src/VectorRandomFeature.jl index 6abcc83dd..1e029c5ea 100644 --- a/src/VectorRandomFeature.jl +++ b/src/VectorRandomFeature.jl @@ -151,14 +151,14 @@ Constructs a `VectorRandomFeatureInterface <: MachineLearningTool` interface for - "n_iteration": number of eki iterations - "scheduler": Learning rate Scheduler (a.k.a. EKP timestepper) Default: DataMisfitController - "cov_sample_multiplier": increase for more samples to estimate covariance matrix in optimization (default 10.0, minimum 0.0) - - "tikhonov": tikhonov regularization parameter if > 0 - "inflation": additive inflation ∈ [0,1] with 0 being no inflation - "train_fraction": e.g. 0.8 (default) means 80:20 train - test split - "n_features_opt": fix the number of features for optimization (default `n_features`, as used for prediction) - "multithread": how to multithread. "ensemble" (default) threads across ensemble members "tullio" threads random feature matrix algebra - "accelerator": use EKP accelerators (default is no acceleration) - "verbose" => false, verbose optimizer statements to check convergence, priors and optimal parameters. - - "cov_correction" => "shrinkage", type of conditioning to improve estimated covariance (Ledoit Wolfe 03), also "nice" for (Vishny, Morzfeld et al. 2024) + - "cov_correction" => "nice": type of conditioning to improve estimated covariance. "shrinkage", "shrinkage_corr" (Ledoit Wolfe 03), "nice" for (Vishny, Morzfeld et al. 2024) + - "overfit" => 1.0: if > 1.0 forcibly overfit/under-regularize the optimizer cost, (vice versa for < 1.0). - "n_cross_val_sets" => 2, train fraction creates (default 5) train-test data subsets, then use 'n_cross_val_sets' of these stacked in the loss function. If set to 0, train=test on the full data provided ignoring "train_fraction". """ @@ -194,7 +194,6 @@ function VectorRandomFeatureInterface( "n_iteration" => 10, # number of eki iterations "scheduler" => EKP.DataMisfitController(terminate_at = 1000), # Adaptive timestepping "cov_sample_multiplier" => 10.0, # multiplier for samples to estimate covariance in optimization scheme - "tikhonov" => 0, # tikhonov regularization parameter if >0 "inflation" => 1e-4, # additive inflation ∈ [0,1] with 0 being no inflation "train_fraction" => 0.8, # 80:20 train - test split "n_features_opt" => n_features, # number of features for the optimization @@ -202,8 +201,9 @@ function VectorRandomFeatureInterface( "verbose" => false, # verbose optimizer statements "localization" => EKP.Localizers.NoLocalization(), # localization / sample error correction for small ensembles "accelerator" => EKP.NesterovAccelerator(), # acceleration with momentum - "cov_correction" => "shrinkage", # type of conditioning to improve estimated covariance + "cov_correction" => "nice", # type of conditioning to improve estimated covariance "n_cross_val_sets" => 2, # if set to 0, removes data split. i.e takes train & test to be the same data set + "overfit" => 1.0, # if >1 this forcibly overfits to the data ) if !isnothing(optimizer_options) @@ -218,7 +218,9 @@ function VectorRandomFeatureInterface( opt_tmp[key] = optimizer_opts[key] end end - @info("hyperparameter optimization with EKI configured with $opt_tmp") + if optimizer_opts["verbose"] + @info("hyperparameter optimization with EKI configured with $opt_tmp") + end return VectorRandomFeatureInterface{S, RNG, KSType}( rfms, @@ -243,21 +245,23 @@ function hyperparameter_distribution_from_flat( input_dim::Int, output_dim::Int, kernel_structure::SK, + prior_in_scale, ) where {VV <: AbstractVector, SK <: SeparableKernel} - M = zeros(input_dim, output_dim) # n x p mean U, V = hyperparameters_from_flat(x, input_dim, output_dim, kernel_structure) - - if !isposdef(U) - println("U not posdef - correcting") - U = posdef_correct(U) + Uscaled = Diagonal(vec(prior_in_scale)) * U * Diagonal(vec(prior_in_scale)) # scale U only + Uscaled = 0.5 * (Uscaled + Uscaled') + V = 0.5 * (V + V') + if !isposdef(Uscaled) + println("U not posdef - correcting. Consider increasing argument eps in the chosen kernel structure") + Uscaled = posdef_correct(Uscaled) end if !isposdef(V) - println("V not posdef - correcting") + println("V not posdef - correcting. Consider increasing argument eps in the chosen kernel structure") V = posdef_correct(V) end - dist = MatrixNormal(M, U, V) + dist = MatrixNormal(zeros(input_dim, output_dim), Uscaled, V) pd = ParameterDistribution( Dict( "distribution" => Parameterized(dist), @@ -273,14 +277,18 @@ function hyperparameter_distribution_from_flat( input_dim::Int, output_dim::Int, kernel_structure::NK, + prior_in_scale, ) where {VV <: AbstractVector, NK <: NonseparableKernel} C = hyperparameters_from_flat(x, input_dim, output_dim, kernel_structure) - if !isposdef(C) - println("C not posdef - correcting") - C = posdef_correct(C) + scale_vec = reduce(vcat, repeat(vec(prior_in_scale), output_dim)) # make repeating vector ("out_dim" times) + Cscaled = Diagonal(vec(scale_vec)) * C * Diagonal(vec(scale_vec)) + Cscaled = 0.5 * (Cscaled + Cscaled') + if !isposdef(Cscaled) + println("C not posdef - correcting, consider increasing argument eps in the chosen kernel structure") + Cscaled = posdef_correct(Cscaled) end - dist = MvNormal(zeros(input_dim * output_dim), C) + dist = MvNormal(zeros(input_dim * output_dim), Cscaled) pd = ParameterDistribution( Dict( "distribution" => Parameterized(dist), @@ -306,7 +314,9 @@ function RFM_from_hyperparameters( batch_sizes::Union{Dict{S, Int}, Nothing}, input_dim::Int, output_dim::Int, - multithread_type::MT; + multithread_type::MT, + prior_in_scale, + prior_out_scale, ) where { S <: AbstractString, RNG <: AbstractRNG, @@ -315,15 +325,15 @@ function RFM_from_hyperparameters( MT <: MultithreadType, } - xi_hp = l[1:(end - 1)] - sigma_hp = l[end] + xi_hp = isa(l, AbstractVecOrMat) ? l[:] : [l] + prior_out_scale_scalar = maximum(prior_out_scale) # most conservative scaling kernel_structure = get_kernel_structure(vrfi) - pd = hyperparameter_distribution_from_flat(xi_hp, input_dim, output_dim, kernel_structure) + pd = hyperparameter_distribution_from_flat(xi_hp, input_dim, output_dim, kernel_structure, prior_in_scale) feature_sampler = RF.Samplers.FeatureSampler(pd, output_dim, rng = rng) - feature_parameters = Dict("sigma" => sigma_hp) + feature_parameters = Dict("sigma" => prior_out_scale_scalar * sqrt(2)) vff = RF.Features.VectorFourierFeature( n_features, output_dim, @@ -380,6 +390,7 @@ function build_models!( batch_sizes = get_batch_sizes(vrfi) decomp_type = get_feature_decomposition(vrfi) optimizer_options = get_optimizer_options(vrfi) + opt_verbose_flag = optimizer_options["verbose"] optimizer = get_optimizer(vrfi) multithread = optimizer_options["multithread"] if multithread == "ensemble" @@ -393,7 +404,7 @@ function build_models!( ), ) end - prior = optimizer_options["prior"] + prior = build_default_prior(input_dim, output_dim, kernel_structure) rng = get_rng(vrfi) # where prior space has changed we need to rebuild the priors [TODO just build them here in the first place?] @@ -405,6 +416,10 @@ function build_models!( prior = build_default_prior(input_dim, output_dim, kernel_structure) end + # scale up the prior so that default priors are always "reasonable" + prior_in_scale = 1.0 ./ std(input_values, dims = 2) + prior_out_scale = std(output_values, dims = 2) + # Optimize feature cholesky factors with EKP # [1.] Split data into test/train (e.g. 80/20) n_features_opt = optimizer_options["n_features_opt"] @@ -460,6 +475,8 @@ function build_models!( μ_hp = transform_unconstrained_to_constrained(prior, mean(prior)) cov_sample_multiplier = optimizer_options["cov_sample_multiplier"] cov_correction = optimizer_options["cov_correction"] + overfit = max(optimizer_options["overfit"], 1e-4) + if nameof(typeof(kernel_structure)) == :SeparableKernel if nameof(typeof(get_output_cov_structure(kernel_structure))) == :DiagonalFactor n_cov_samples_min = n_test + 2 # diagonal case @@ -471,7 +488,6 @@ function build_models!( end n_cov_samples = Int(floor(n_cov_samples_min * max(cov_sample_multiplier, 0.0))) observation_vec = [] - tikhonov_opt_val = optimizer_options["tikhonov"] for cv_idx in 1:n_cross_val_sets internal_Γ, approx_σ2 = estimate_mean_and_coeffnorm_covariance( vrfi, @@ -486,40 +502,20 @@ function build_models!( n_cov_samples, decomp_type, multithread_type, + prior_in_scale, + prior_out_scale, cov_correction = cov_correction, + verbose = opt_verbose_flag, ) - if tikhonov_opt_val == 0 - # Build the covariance - Γ = internal_Γ - Γ[1:(n_test * output_dim), 1:(n_test * output_dim)] += - isa(regularization, UniformScaling) ? regularization : kron(I(n_test), regularization) # + approx_σ2 - Γ[(n_test * output_dim + 1):end, (n_test * output_dim + 1):end] += I - data = vcat(reshape(get_outputs(input_output_pairs)[:, test_idx[cv_idx]], :, 1), 0.0, 0.0) #flatten data - - elseif tikhonov_opt_val > 0 - # augment the state to add tikhonov - outsize = size(internal_Γ, 1) - Γ = zeros(outsize + n_hp, outsize + n_hp) - Γ[1:outsize, 1:outsize] = internal_Γ - Γ[1:(n_test * output_dim), 1:(n_test * output_dim)] += kron(I(n_test), regularization) # block diag regularization - Γ[(n_test * output_dim + 1):outsize, (n_test * output_dim + 1):outsize] += I - - Γ[(outsize + 1):end, (outsize + 1):end] = tikhonov_opt_val .* cov(prior) - - data = vcat( - reshape(get_outputs(input_output_pairs)[:, test_idx[cv_idx]], :, 1), - 0.0, - 0.0, - zeros(size(Γ, 1) - outsize, 1), - ) #flatten data with additional zeros - else - throw( - ArgumentError( - "Tikhonov parameter must be non-negative, instead received tikhonov_opt_val=$tikhonov_opt_val", - ), - ) - end + # Build the covariance + Γ = deepcopy(internal_Γ) + Γ[1:(n_test * output_dim), 1:(n_test * output_dim)] += + isa(regularization, UniformScaling) ? regularization : kron(I(n_test), regularization) # + approx_σ2 + Γ[1:(n_test * output_dim), 1:(n_test * output_dim)] /= overfit^2 + Γ[(n_test * output_dim + 1):end, (n_test * output_dim + 1):end] += I + data = vcat(reshape(get_outputs(input_output_pairs)[:, test_idx[cv_idx]], :, 1), 0.0, 0.0) #flatten data + if !isposdef(Γ) Γ = posdef_correct(Γ) end @@ -538,11 +534,14 @@ function build_models!( initial_params = construct_initial_ensemble(rng, prior, n_ensemble) + # bug with scalar mean o/w + prior_mean = isa(mean(prior), AbstractVector) ? mean(prior) : [mean(prior)] + prior_cov = cov(prior) ekiobj = EKP.EnsembleKalmanProcess( initial_params, observation, - Inversion(), + TransformInversion(prior_mean, prior_cov), scheduler = scheduler, rng = rng, verbose = opt_verbose_flag, @@ -556,11 +555,7 @@ function build_models!( #get parameters: lvec = get_ϕ_final(prior, ekiobj) - if tikhonov_opt_val > 0 - g_ens = zeros(n_cross_val_sets * (output_dim * n_test + input_dim + 2), n_ensemble) - else - g_ens = zeros(n_cross_val_sets * (output_dim * n_test + 2), n_ensemble) - end + g_ens = zeros(n_cross_val_sets * (output_dim * n_test + 2), n_ensemble) for cv_idx in 1:n_cross_val_sets g_ens_tmp, _ = calculate_ensemble_mean_and_coeffnorm( vrfi, @@ -574,20 +569,13 @@ function build_models!( input_output_pairs, decomp_type, multithread_type, + prior_in_scale, + prior_out_scale, + verbose = opt_verbose_flag, ) - if tikhonov_opt_val > 0 - # augment with the computational parameters (u not ϕ) - uvecormat = get_u_final(ekiobj) - if isa(uvecormat, AbstractVector) - umat = reshape(uvecormat, 1, :) - else - umat = uvecormat - end + indices = ((cv_idx - 1) * (output_dim * n_test + 2) + 1):(cv_idx * (output_dim * n_test + 2)) - g_ens_tmp = vcat(g_ens_tmp, umat) - end - - g_ens[((cv_idx - 1) * (output_dim * n_test + 2) + 1):(cv_idx * (output_dim * n_test + 2)), :] = g_ens_tmp + g_ens[indices, :] = g_ens_tmp end @@ -641,10 +629,11 @@ function build_models!( input_dim, output_dim, multithread_type, + prior_in_scale, + prior_out_scale, ) fitted_features_tmp = fit(rfm_tmp, input_output_pairs, decomposition_type = decomp_type) - push!(rfms, rfm_tmp) push!(fitted_features, fitted_features_tmp) push!(get_regularization(vrfi), regularization) diff --git a/test/RandomFeature/runtests.jl b/test/RandomFeature/runtests.jl index 90feb3a38..3778e8c25 100644 --- a/test/RandomFeature/runtests.jl +++ b/test/RandomFeature/runtests.jl @@ -45,7 +45,7 @@ rng = Random.MersenneTwister(seed) @test calculate_n_hyperparameters(d, odf) == 0 @test calculate_n_hyperparameters(d, df) == d @test calculate_n_hyperparameters(d, cf) == Int(d * (d + 1) / 2) - @test calculate_n_hyperparameters(d, lrf) == Int(r + d * r) + @test calculate_n_hyperparameters(d, lrf) == Int(d + d * r) @test calculate_n_hyperparameters(d, hlrf) == Int(r * (r + 1) / 2 + d * r) # hyperparameters_from_flat - only check shape @@ -64,6 +64,22 @@ rng = Random.MersenneTwister(seed) @test ndims(prior) == n_hp end + # API from string: + d = 12 + cov_strings_and_checks = [ + ("onedim", OneDimFactor()), + ("diagonal", DiagonalFactor()), + ("cholesky", CholeskyFactor()), + ("lowrank", LowRankFactor(Int(ceil(sqrt(d))))), + ("hierlowrank", HierarchicalLowRankFactor(Int(ceil(sqrt(d))))), + ] + + for (cs, cc) in cov_strings_and_checks + @test cov_structure_from_string(cs, d) == cc + end + @test cov_structure_from_string(OneDimFactor()) == OneDimFactor() + @test_throws ArgumentError cov_structure_from_string("bad-string", d) + # [2. ] Kernel Structures d = 6 p = 3 @@ -79,12 +95,12 @@ rng = Random.MersenneTwister(seed) # calculate_n_hyperparameters @test calculate_n_hyperparameters(d, p, k_sep) == - calculate_n_hyperparameters(d, c_in) + calculate_n_hyperparameters(p, c_out) + 1 - @test calculate_n_hyperparameters(d, p, k_nonsep) == calculate_n_hyperparameters(d * p, c_in) + 1 + calculate_n_hyperparameters(d, c_in) + calculate_n_hyperparameters(p, c_out) + @test calculate_n_hyperparameters(d, p, k_nonsep) == calculate_n_hyperparameters(d * p, c_in) k_sep1d = SeparableKernel(odf, odf) k_nonsep1d = NonseparableKernel(odf) - @test calculate_n_hyperparameters(1, 1, k_sep1d) == 2 - @test calculate_n_hyperparameters(1, 1, k_nonsep1d) == 2 + @test calculate_n_hyperparameters(1, 1, k_sep1d) == 1 + @test calculate_n_hyperparameters(1, 1, k_nonsep1d) == 1 # hyper_parameters_from_flat: not applied with scaling hyperparameter x = ones(calculate_n_hyperparameters(d, c_in) + calculate_n_hyperparameters(p, c_out)) @@ -102,28 +118,29 @@ rng = Random.MersenneTwister(seed) # build_default_prior @test ndims(build_default_prior(d, p, k_sep)) == ndims(build_default_prior("input", calculate_n_hyperparameters(d, c_in), c_in)) + - ndims(build_default_prior("output", calculate_n_hyperparameters(p, c_out), c_out)) + - 1 + ndims(build_default_prior("output", calculate_n_hyperparameters(p, c_out), c_out)) @test ndims(build_default_prior(d, p, k_nonsep)) == - ndims(build_default_prior("full", calculate_n_hyperparameters(d * p, c_in), c_in)) + 1 + ndims(build_default_prior("full", calculate_n_hyperparameters(d * p, c_in), c_in)) - @test ndims(build_default_prior(1, 1, k_sep1d)) == 2 - @test ndims(build_default_prior(1, 1, k_nonsep1d)) == 2 + @test ndims(build_default_prior(1, 1, k_sep1d)) == 1 + @test ndims(build_default_prior(1, 1, k_nonsep1d)) == 1 # test shrinkage utility - samples = rand(MvNormal(zeros(100), I), 20) + samples = rand(MvNormal(zeros(100), I), 50) # normal condition number should be huge around 10^18 # shrinkage cov will have condition number around 1 and be close to I good_cov = shrinkage_cov(samples) - @test (cond(good_cov) < 10) && ((good_cov[1] < 1.5) && (good_cov[1] > 0.5)) + @test (cond(good_cov) < 10) && ((good_cov[1] < 5.0) && (good_cov[1] > 0.2)) + + good_cov = shrinkage_cov(samples, cov_or_corr = "corr", verbose = true) + @test (cond(good_cov) < 10) && ((good_cov[1] < 5.0) && (good_cov[1] > 0.2)) # test NICE utility - samples = rand(MvNormal(zeros(100), I), 20) # normal condition number should be huge around 10^18 # nice cov will have improved conditioning, does not perform as well at this task as shrinking so has looser bounds - good_cov = nice_cov(samples) - @test (cond(good_cov) < 100) && ((good_cov[1] < 2.0) && (good_cov[1] > 0.2)) + good_cov = nice_cov(samples, verbose = true) + @test (cond(good_cov) < 100) && ((good_cov[1] < 5.0) && (good_cov[1] > 0.2)) end @@ -152,7 +169,7 @@ rng = Random.MersenneTwister(seed) "multithread" => "ensemble", "accelerator" => NesterovAccelerator(), "verbose" => false, - "cov_correction" => "shrinkage", + "cov_correction" => "nice", "n_cross_val_sets" => 2, ) @@ -210,15 +227,15 @@ rng = Random.MersenneTwister(seed) "scheduler" => DataMisfitController(terminate_at = 1000), "cov_sample_multiplier" => 10.0, "n_features_opt" => n_features, - "tikhonov" => 0, "inflation" => 1e-4, "train_fraction" => 0.8, "multithread" => "ensemble", "accelerator" => NesterovAccelerator(), "verbose" => false, "localization" => EnsembleKalmanProcesses.Localizers.NoLocalization(), - "cov_correction" => "shrinkage", + "cov_correction" => "nice", "n_cross_val_sets" => 2, + "overfit" => 1.0, ) #build interfaces @@ -286,7 +303,6 @@ rng = Random.MersenneTwister(seed) eps = 1e-8 # more reg needed here for some reason... vector_ks = SeparableKernel(DiagonalFactor(eps), CholeskyFactor()) # Diagonalize input (ARD-type kernel) # Scalar RF options to mimic squared-exp ARD kernel - n_features = 100 srfi = ScalarRandomFeatureInterface( n_features, input_dim, @@ -389,25 +405,32 @@ rng = Random.MersenneTwister(seed) # RF parameters - n_features = 150 - # Test a few options for RF + # Test a few options branches for RF # 1) scalar + diag in # 2) scalar - # 3) vector + diag out - # 4) vector - # 5) vector nonseparable - eps = 1e-8 - r = 1 + # 3) vector + diag out, correct cov by shrinkage (cov) + # 4) vector , correct cov by shrinkage (corr) + # 5) vector nonseparable , default correction with "nice" + eps = 1e-6 + r_sep = 1 + r_nonsep = 2 scalar_diagin_ks = SeparableKernel(DiagonalFactor(eps), OneDimFactor()) scalar_ks = SeparableKernel(CholeskyFactor(eps), OneDimFactor()) vector_diagout_ks = SeparableKernel(CholeskyFactor(eps), DiagonalFactor(eps)) - vector_ks = SeparableKernel(LowRankFactor(r, eps), HierarchicalLowRankFactor(r, eps)) - vector_nonsep_ks = NonseparableKernel(LowRankFactor(r + 1, eps)) + vector_ks = SeparableKernel(HierarchicalLowRankFactor(r_sep, eps), LowRankFactor(r_sep, eps)) + vector_nonsep_ks = NonseparableKernel(LowRankFactor(r_nonsep, eps)) + n_features = 100 srfi_diagin = ScalarRandomFeatureInterface(n_features, input_dim, kernel_structure = scalar_diagin_ks, rng = rng) - srfi = ScalarRandomFeatureInterface(n_features, input_dim, kernel_structure = scalar_ks, rng = rng) + srfi = ScalarRandomFeatureInterface( + n_features, + input_dim, + kernel_structure = scalar_ks, + rng = rng, + optimizer_options = Dict("verbose" => true), + ) vrfi_diagout = VectorRandomFeatureInterface( n_features, @@ -415,9 +438,17 @@ rng = Random.MersenneTwister(seed) output_dim, kernel_structure = vector_diagout_ks, rng = rng, + optimizer_options = Dict("cov_correction" => "shrinkage_corr"), ) - vrfi = VectorRandomFeatureInterface(n_features, input_dim, output_dim, kernel_structure = vector_ks, rng = rng) + vrfi = VectorRandomFeatureInterface( + n_features, + input_dim, + output_dim, + kernel_structure = vector_ks, + rng = rng, + optimizer_options = Dict("verbose" => true, "cov_correction" => "shrinkage"), + ) vrfi_nonsep = VectorRandomFeatureInterface( n_features, @@ -425,6 +456,7 @@ rng = Random.MersenneTwister(seed) output_dim, kernel_structure = vector_nonsep_ks, rng = rng, + optimizer_options = Dict("verbose" => true), ) # build emulators @@ -434,9 +466,9 @@ rng = Random.MersenneTwister(seed) em_srfi_svd = Emulator(srfi, iopairs; encoder_kwargs = (; obs_noise_cov = Σ)) em_vrfi_svd_diagout = Emulator(vrfi_diagout, iopairs; encoder_kwargs = (; obs_noise_cov = Σ)) - em_vrfi_svd = Emulator(vrfi, iopairs; encoder_kwargs = (; obs_noise_cov = Σ)) + em_vrfi_svd = Emulator(deepcopy(vrfi), iopairs; encoder_kwargs = (; obs_noise_cov = Σ)) - em_vrfi = Emulator(vrfi, iopairs; encoder_schedule = [], encoder_kwargs = (; obs_noise_cov = Σ)) + em_vrfi = Emulator(deepcopy(vrfi), iopairs; encoder_schedule = [], encoder_kwargs = (; obs_noise_cov = Σ)) em_vrfi_nonsep = Emulator(vrfi_nonsep, iopairs; encoder_schedule = [], encoder_kwargs = (; obs_noise_cov = Σ)) #TODO truncated SVD option for vector (involves resizing prior) @@ -460,16 +492,15 @@ rng = Random.MersenneTwister(seed) @test isapprox.(norm(μ_ss - outx), 0, atol = tol_μ) @test isapprox.(norm(μ_vsd - outx), 0, atol = tol_μ) @test isapprox.(norm(μ_vs - outx), 0, atol = tol_μ) - @test isapprox.(norm(μ_v - outx), 0, atol = 10 * tol_μ) # approximate option so likely less good approx - @test isapprox.(norm(μ_vns - outx), 0, atol = 10 * tol_μ) # approximate option so likely less good approx - + @test isapprox.(norm(μ_v - outx), 0, atol = 2 * tol_μ) # approximate option so likely less good approx + @test isapprox.(norm(μ_vns - outx), 0, atol = 2 * tol_μ) # approximate option so likely less good approx # An example with the other threading option vrfi_tul = VectorRandomFeatureInterface( n_features, input_dim, output_dim, kernel_structure = vector_ks, - optimizer_options = Dict("train_fraction" => 0.7, "multithread" => "tullio"), + optimizer_options = Dict("train_fraction" => 0.8, "multithread" => "tullio"), rng = rng, ) em_vrfi_svd_tul = Emulator(vrfi_tul, iopairs; encoder_kwargs = (; obs_noise_cov = Σ))