diff --git a/.gitignore b/.gitignore index 83b96794c..60388f1eb 100644 --- a/.gitignore +++ b/.gitignore @@ -15,4 +15,5 @@ docs/site/ *.vscode* *.DS_Store *.jld2 +*~ output* diff --git a/Project.toml b/Project.toml index 2812a22b1..8a72abcef 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "CalibrateEmulateSample" uuid = "95e48a1f-0bec-4818-9538-3db4340308e3" -authors = ["CLIMA contributors "] version = "0.7.0" +authors = ["CLIMA contributors "] [deps] AbstractGPs = "99985d1d-32ba-4be9-9821-2ec096f28918" @@ -16,6 +16,8 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" GaussianProcesses = "891a1506-143c-57d2-908e-e1f8e92e6de9" KernelFunctions = "ec8451be-7e33-11e9-00cf-bbf324bd1392" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e" +LowRankApprox = "898213cb-b102-5a47-900c-97e73b919f73" MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" @@ -28,6 +30,7 @@ ScikitLearn = "3646fa90-6ef7-5e7e-9f22-8aca16db6324" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e" [compat] AbstractGPs = "0.5.21" @@ -41,7 +44,10 @@ EnsembleKalmanProcesses = "2" ForwardDiff = "1" GaussianProcesses = "0.12" KernelFunctions = "0.10.64" +LinearMaps = "3.11.4" +LowRankApprox = "0.5.5" MCMCChains = "7" +Plots = "1.41.1" Printf = "1" ProgressBars = "1" PyCall = "1.93" @@ -52,6 +58,7 @@ ScikitLearn = "0.6, 0.7" StableRNGs = "1" Statistics = "1" StatsBase = "0.33, 0.34" +TSVD = "0.4.4" julia = "1.6, 1.7, 1.8, 1.9, 1.10, 1.11" [extras] diff --git a/docs/Manifest.toml b/docs/Manifest.toml index 210f1e18e..188bad891 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -5,9 +5,9 @@ manifest_format = "2.0" project_hash = "e4d965798d55f5903483b71f545533e62ea32fdf" [[deps.ADTypes]] -git-tree-sha1 = "be7ae030256b8ef14a441726c4c37766b90b93a3" +git-tree-sha1 = "8be2ae325471fc20b11c27bb34b518541d07dd3a" uuid = "47edcb42-4c32-4615-8424-f2b9edc5f35b" -version = "1.15.0" +version = "1.19.0" [deps.ADTypes.extensions] ADTypesChainRulesCoreExt = "ChainRulesCore" @@ -60,9 +60,9 @@ version = "0.4.5" [[deps.Adapt]] deps = ["LinearAlgebra", "Requires"] -git-tree-sha1 = "f7817e2e585aa6d924fd714df1e2a84be7896c60" +git-tree-sha1 = "7e35fca2bdfba44d797c53dfe63a51fabf39bfc0" uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -version = "4.3.0" +version = "4.4.0" weakdeps = ["SparseArrays", "StaticArrays"] [deps.Adapt.extensions] @@ -115,18 +115,19 @@ version = "3.5.1+1" [[deps.ArrayInterface]] deps = ["Adapt", "LinearAlgebra"] -git-tree-sha1 = "9606d7832795cbef89e06a550475be300364a8aa" +git-tree-sha1 = "d81ae5489e13bc03567d4fbbb06c546a5e53c857" uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -version = "7.19.0" +version = "7.22.0" [deps.ArrayInterface.extensions] ArrayInterfaceBandedMatricesExt = "BandedMatrices" ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices" ArrayInterfaceCUDAExt = "CUDA" - ArrayInterfaceCUDSSExt = "CUDSS" + ArrayInterfaceCUDSSExt = ["CUDSS", "CUDA"] ArrayInterfaceChainRulesCoreExt = "ChainRulesCore" ArrayInterfaceChainRulesExt = "ChainRules" ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore" + ArrayInterfaceMetalExt = "Metal" ArrayInterfaceReverseDiffExt = "ReverseDiff" ArrayInterfaceSparseArraysExt = "SparseArrays" ArrayInterfaceStaticArraysCoreExt = "StaticArraysCore" @@ -140,6 +141,7 @@ version = "7.19.0" ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" + Metal = "dde4c033-4e86-420c-a63e-0dd931031962" ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" @@ -157,9 +159,9 @@ version = "1.1.0" [[deps.AxisArrays]] deps = ["Dates", "IntervalSets", "IterTools", "RangeArrays"] -git-tree-sha1 = "16351be62963a67ac4083f748fdb3cca58bfd52f" +git-tree-sha1 = "4126b08903b777c88edf1754288144a0492c05ad" uuid = "39de3d68-74b9-583c-8d2d-e117c070f3a9" -version = "0.4.7" +version = "0.4.8" [[deps.BangBang]] deps = ["Compat", "ConstructionBase", "InitialValues", "LinearAlgebra", "Requires", "Setfield", "Tables"] @@ -192,9 +194,9 @@ version = "0.1.1" [[deps.BenchmarkTools]] deps = ["Compat", "JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"] -git-tree-sha1 = "e38fbc49a620f5d0b660d7f543db1009fe0f8336" +git-tree-sha1 = "7fecfb1123b8d0232218e2da0c213004ff15358d" uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" -version = "1.6.0" +version = "1.6.3" [[deps.BitTwiddlingConvenienceFunctions]] deps = ["Static"] @@ -209,10 +211,10 @@ uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0" version = "1.0.9+0" [[deps.CPUSummary]] -deps = ["CpuId", "IfElse", "PrecompileTools", "Static"] -git-tree-sha1 = "5a97e67919535d6841172016c9530fd69494e5ec" +deps = ["CpuId", "IfElse", "PrecompileTools", "Preferences", "Static"] +git-tree-sha1 = "f3a21d7fc84ba618a779d1ed2fcca2e682865bab" uuid = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9" -version = "0.2.6" +version = "0.2.7" [[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"] @@ -221,16 +223,16 @@ uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a" 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"] +deps = ["AbstractGPs", "AbstractMCMC", "AdvancedMH", "ChunkSplitters", "Conda", "Distributions", "DocStringExtensions", "EnsembleKalmanProcesses", "ForwardDiff", "GaussianProcesses", "KernelFunctions", "LinearAlgebra", "LinearMaps", "LowRankApprox", "MCMCChains", "Pkg", "Printf", "ProgressBars", "PyCall", "Random", "RandomFeatures", "ReverseDiff", "ScikitLearn", "StableRNGs", "Statistics", "StatsBase", "TSVD"] path = ".." uuid = "95e48a1f-0bec-4818-9538-3db4340308e3" version = "0.7.0" [[deps.ChainRulesCore]] deps = ["Compat", "LinearAlgebra"] -git-tree-sha1 = "06ee8d1aa558d2833aa799f6f0b31b30cada405f" +git-tree-sha1 = "e4c6a16e77171a5f5e25e9646617ab1c276c5607" uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "1.25.2" +version = "1.26.0" weakdeps = ["SparseArrays"] [deps.ChainRulesCore.extensions] @@ -272,9 +274,9 @@ version = "1.0.0" [[deps.Compat]] deps = ["TOML", "UUIDs"] -git-tree-sha1 = "3a3dfb30697e96a440e4149c8c51bf32f818c0f3" +git-tree-sha1 = "9d8a54ce4b17aa5bdce0ea5c34bc5e7c340d16ad" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "4.17.0" +version = "4.18.1" weakdeps = ["Dates", "LinearAlgebra"] [deps.Compat.extensions] @@ -296,9 +298,9 @@ weakdeps = ["InverseFunctions"] [[deps.Conda]] deps = ["Downloads", "JSON", "VersionParsing"] -git-tree-sha1 = "b19db3927f0db4151cb86d073689f2428e524576" +git-tree-sha1 = "8f06b0cfa4c514c7b9546756dbae91fcfbc92dc9" uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d" -version = "1.10.2" +version = "1.10.3" [[deps.ConsoleProgressMonitor]] deps = ["Logging", "ProgressMeter"] @@ -341,9 +343,9 @@ version = "1.16.0" [[deps.DataFrames]] deps = ["Compat", "DataAPI", "DataStructures", "Future", "InlineStrings", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrecompileTools", "PrettyTables", "Printf", "Random", "Reexport", "SentinelArrays", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"] -git-tree-sha1 = "fb61b4812c49343d7ef0b533ba982c46021938a6" +git-tree-sha1 = "d8928e9169ff76c6281f39a659f9bca3a573f24c" uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" -version = "1.7.0" +version = "1.8.1" [[deps.DataStructures]] deps = ["Compat", "InteractiveUtils", "OrderedCollections"] @@ -380,9 +382,9 @@ version = "1.15.1" [[deps.DifferentiationInterface]] deps = ["ADTypes", "LinearAlgebra"] -git-tree-sha1 = "210933c93f39f832d92f9efbbe69a49c453db36d" +git-tree-sha1 = "80bd15222b3e8d0bc70d921d2201aa0084810ce5" uuid = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63" -version = "0.7.1" +version = "0.7.12" [deps.DifferentiationInterface.extensions] DifferentiationInterfaceChainRulesCoreExt = "ChainRulesCore" @@ -446,9 +448,9 @@ version = "1.11.0" [[deps.Distributions]] deps = ["AliasTables", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns"] -git-tree-sha1 = "3e6d038b77f22791b8e3472b7c633acea1ecac06" +git-tree-sha1 = "3bc002af51045ca3b47d2e1787d6ce02e68b943a" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" -version = "0.25.120" +version = "0.25.122" [deps.Distributions.extensions] DistributionsChainRulesCoreExt = "ChainRulesCore" @@ -467,9 +469,9 @@ version = "0.9.5" [[deps.Documenter]] 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" +git-tree-sha1 = "70c521ca3a23c576e12655d15963977c9766c26b" uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -version = "1.13.0" +version = "1.16.0" [[deps.Downloads]] deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] @@ -484,15 +486,21 @@ version = "1.2.12" [[deps.ElasticPDMats]] deps = ["LinearAlgebra", "MacroTools", "PDMats"] -git-tree-sha1 = "03ec11d0151e8a772b396aecd663e1c76fc8edcf" +git-tree-sha1 = "feb325f9ca96dce03ae03347ff23854205f5be3d" uuid = "2904ab23-551e-5aed-883f-487f97af5226" -version = "0.2.3" +version = "0.2.4" [[deps.EnsembleKalmanProcesses]] deps = ["Convex", "Distributions", "DocStringExtensions", "FFMPEG", "GaussianRandomFields", "Interpolations", "LinearAlgebra", "MathOptInterface", "Optim", "QuadGK", "Random", "RecipesBase", "SCS", "SparseArrays", "Statistics", "StatsBase", "TOML", "TSVD"] -git-tree-sha1 = "dd6b4b258beeef8db8024e14fea95d34305e5bb0" +git-tree-sha1 = "52faab1c57dac15c161f60f0c1d94ebb110ea3c8" uuid = "aa8a2aa5-91d8-4396-bcef-d4f2ec43552d" -version = "2.4.0" +version = "2.6.0" + + [deps.EnsembleKalmanProcesses.extensions] + EnsembleKalmanProcessesMakieExt = "Makie" + + [deps.EnsembleKalmanProcesses.weakdeps] + Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" [[deps.EnumX]] git-tree-sha1 = "bddad79635af6aec424f53ed8aad5d7555dc6f00" @@ -501,27 +509,27 @@ version = "1.0.5" [[deps.Expat_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "d55dffd9ae73ff72f1c0482454dcf2ec6c6c4a63" +git-tree-sha1 = "27af30de8b5445644e8ffe3bcb0d72049c089cf1" uuid = "2e619515-83b5-522b-bb60-26c02a35a201" -version = "2.6.5+0" +version = "2.7.3+0" [[deps.FFMPEG]] deps = ["FFMPEG_jll"] -git-tree-sha1 = "53ebe7511fa11d33bec688a9178fac4e49eeee00" +git-tree-sha1 = "95ecf07c2eea562b5adbd0696af6db62c0f52560" uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a" -version = "0.4.2" +version = "0.4.5" [[deps.FFMPEG_jll]] deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "PCRE2_jll", "Zlib_jll", "libaom_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"] -git-tree-sha1 = "466d45dc38e15794ec7d5d63ec03d776a9aff36e" +git-tree-sha1 = "ccc81ba5e42497f4e76553a5545665eed577a663" uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5" -version = "4.4.4+1" +version = "8.0.0+0" [[deps.FFTW]] -deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"] -git-tree-sha1 = "797762812ed063b9b94f6cc7742bc8883bb5e69e" +deps = ["AbstractFFTs", "FFTW_jll", "Libdl", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"] +git-tree-sha1 = "97f08406df914023af55ade2f843c39e99c5d969" uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" -version = "1.9.0" +version = "1.10.0" [[deps.FFTW_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -541,9 +549,9 @@ version = "1.11.0" [[deps.FillArrays]] deps = ["LinearAlgebra"] -git-tree-sha1 = "6a70198746448456524cb442b8af316927ff3e1a" +git-tree-sha1 = "5bfcd42851cf2f1b303f51525a54dc5e98d408a3" uuid = "1a297f60-69ca-5386-bcde-b61e274b549b" -version = "1.13.0" +version = "1.15.0" weakdeps = ["PDMats", "SparseArrays", "Statistics"] [deps.FillArrays.extensions] @@ -553,9 +561,9 @@ weakdeps = ["PDMats", "SparseArrays", "Statistics"] [[deps.FiniteDiff]] deps = ["ArrayInterface", "LinearAlgebra", "Setfield"] -git-tree-sha1 = "f089ab1f834470c525562030c8cfde4025d5e915" +git-tree-sha1 = "9340ca07ca27093ff68418b7558ca37b05f8aeb1" uuid = "6a86dc24-6348-571c-b903-95158fe2bd41" -version = "2.27.0" +version = "2.29.0" [deps.FiniteDiff.extensions] FiniteDiffBandedMatricesExt = "BandedMatrices" @@ -571,9 +579,9 @@ version = "2.27.0" [[deps.Fontconfig_jll]] deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Zlib_jll"] -git-tree-sha1 = "301b5d5d731a0654825f1f2e906990f7141a106b" +git-tree-sha1 = "f85dac9a96a01087df6e3a749840015a0ca3817d" uuid = "a3f928ae-7b40-5064-980b-68af3947d34b" -version = "2.16.0+0" +version = "2.17.1+0" [[deps.Formatting]] deps = ["Logging", "Printf"] @@ -583,9 +591,9 @@ version = "0.4.3" [[deps.ForwardDiff]] deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"] -git-tree-sha1 = "a2df1b776752e3f344e5116c06d75a10436ab853" +git-tree-sha1 = "afb7c51ac63e40708a3071f80f5e84a752299d4f" uuid = "f6369f11-7733-5829-9624-2563aa707210" -version = "0.10.38" +version = "0.10.39" weakdeps = ["StaticArrays"] [deps.ForwardDiff.extensions] @@ -626,10 +634,10 @@ uuid = "891a1506-143c-57d2-908e-e1f8e92e6de9" version = "0.12.5" [[deps.GaussianRandomFields]] -deps = ["Arpack", "FFTW", "FastGaussQuadrature", "LinearAlgebra", "RecipesBase", "SpecialFunctions", "Statistics"] -git-tree-sha1 = "d9c335f2c06424029b2addf9abf602e0feb2f53e" +deps = ["Arpack", "FFTW", "FastGaussQuadrature", "LinearAlgebra", "Random", "RecipesBase", "SpecialFunctions", "Statistics", "StatsBase"] +git-tree-sha1 = "90f9cad110814a408c66ec7e37cfe61d579af7be" uuid = "e4b2fa32-6e09-5554-b718-106ed5adafe9" -version = "2.1.6" +version = "2.2.7" [[deps.GettextRuntime_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll"] @@ -638,22 +646,28 @@ uuid = "b0724c58-0f36-5564-988d-3bb0596ebc4a" version = "0.22.4+0" [[deps.Git]] -deps = ["Git_jll", "JLLWrappers", "OpenSSH_jll"] -git-tree-sha1 = "2230a9cc32394b11a3b3aa807a382e3bbab1198c" +deps = ["Git_LFS_jll", "Git_jll", "JLLWrappers", "OpenSSH_jll"] +git-tree-sha1 = "824a1890086880696fc908fe12a17bcf61738bd8" uuid = "d7ba0133-e1db-5d97-8f8c-041e4b3a1eb2" -version = "1.4.0" +version = "1.5.0" + +[[deps.Git_LFS_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "bb8471f313ed941f299aa53d32a94ab3bee08844" +uuid = "020c3dae-16b3-5ae5-87b3-4cb189e250b2" +version = "3.7.0+0" [[deps.Git_jll]] deps = ["Artifacts", "Expat_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Libiconv_jll", "OpenSSL_jll", "PCRE2_jll", "Zlib_jll"] -git-tree-sha1 = "b981ed24de5855f20fce5b8cb767c179f93e4268" +git-tree-sha1 = "b6a684587ebe896d9f68ae777f648205940f0f70" uuid = "f8c6e375-362e-5223-8a59-34ff63f689eb" -version = "2.50.0+0" +version = "2.51.3+0" [[deps.Glib_jll]] deps = ["Artifacts", "GettextRuntime_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Zlib_jll"] -git-tree-sha1 = "35fbd0cefb04a516104b8e183ce0df11b70a3f1a" +git-tree-sha1 = "50c11ffab2a3d50192a228c313f05b5b5dc5acb2" uuid = "7746bdde-850d-59dc-9ae8-88ece973131d" -version = "2.84.3+0" +version = "2.86.0+0" [[deps.Graphite2_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -675,9 +689,9 @@ version = "0.1.17" [[deps.IOCapture]] deps = ["Logging", "Random"] -git-tree-sha1 = "b6d6bfdd7ce25b0f9b2f6b3dd56b2673a66c8770" +git-tree-sha1 = "0ee181ec08df7d7c911901ea38baf16f755114dc" uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89" -version = "0.2.5" +version = "1.0.0" [[deps.IfElse]] git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1" @@ -690,9 +704,9 @@ uuid = "22cec73e-a1b8-11e9-2c92-598750a2cf9c" version = "0.3.1" [[deps.InlineStrings]] -git-tree-sha1 = "8594fac023c5ce1ef78260f24d1ad18b4327b420" +git-tree-sha1 = "8f3d257792a522b4601c24a577954b0a8cd7334d" uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48" -version = "1.4.4" +version = "1.4.5" [deps.InlineStrings.extensions] ArrowTypesExt = "ArrowTypes" @@ -704,9 +718,9 @@ version = "1.4.4" [[deps.IntelOpenMP_jll]] deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl"] -git-tree-sha1 = "0f14a5456bdc6b9731a5682f439a672750a09e48" +git-tree-sha1 = "ec1debd61c300961f98064cfb21287613ad7f303" uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0" -version = "2025.0.4+0" +version = "2025.2.0+0" [[deps.InteractiveUtils]] deps = ["Markdown"] @@ -726,12 +740,13 @@ version = "0.15.1" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [[deps.IntervalSets]] -git-tree-sha1 = "5fbb102dcb8b1a858111ae81d56682376130517d" +git-tree-sha1 = "03b4f40b4987baa6a653a21f6f33f902af6255f3" uuid = "8197267c-284f-5f27-9208-e0e47529a953" -version = "0.7.11" -weakdeps = ["Random", "RecipesBase", "Statistics"] +version = "0.7.12" +weakdeps = ["Printf", "Random", "RecipesBase", "Statistics"] [deps.IntervalSets.extensions] + IntervalSetsPrintfExt = "Printf" IntervalSetsRandomExt = "Random" IntervalSetsRecipesBaseExt = "RecipesBase" IntervalSetsStatisticsExt = "Statistics" @@ -768,15 +783,21 @@ version = "1.0.0" [[deps.JLLWrappers]] deps = ["Artifacts", "Preferences"] -git-tree-sha1 = "a007feb38b422fbdab534406aeca1b86823cb4d6" +git-tree-sha1 = "0533e564aae234aff59ab625543145446d8b6ec2" uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" -version = "1.7.0" +version = "1.7.1" [[deps.JSON]] -deps = ["Dates", "Mmap", "Parsers", "Unicode"] -git-tree-sha1 = "31e996f0a15c7b280ba9f76636b3ff9e2ae58c9a" +deps = ["Dates", "Logging", "Parsers", "PrecompileTools", "StructUtils", "UUIDs", "Unicode"] +git-tree-sha1 = "5b6bb73f555bc753a6153deec3717b8904f5551c" uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" -version = "0.21.4" +version = "1.3.0" + + [deps.JSON.extensions] + JSONArrowExt = ["ArrowTypes"] + + [deps.JSON.weakdeps] + ArrowTypes = "31f734f8-188a-4ce0-8406-c8a06bd891cd" [[deps.JSON3]] deps = ["Dates", "Mmap", "Parsers", "PrecompileTools", "StructTypes", "UUIDs"] @@ -792,15 +813,15 @@ version = "1.14.3" [[deps.KernelDensity]] deps = ["Distributions", "DocStringExtensions", "FFTW", "Interpolations", "StatsBase"] -git-tree-sha1 = "7d703202e65efa1369de1279c162b915e245eed1" +git-tree-sha1 = "ba51324b894edaf1df3ab16e2cc6bc3280a2f1a7" uuid = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b" -version = "0.6.9" +version = "0.6.10" [[deps.KernelFunctions]] deps = ["ChainRulesCore", "Compat", "CompositionsBase", "Distances", "FillArrays", "Functors", "IrrationalConstants", "LinearAlgebra", "LogExpFunctions", "Random", "Requires", "SpecialFunctions", "Statistics", "StatsBase", "TensorCore", "Test", "ZygoteRules"] -git-tree-sha1 = "0b8ef8b51580b0d87d0b7a5233bb8ea6d948feb4" +git-tree-sha1 = "0b75b447ee242254ff037402461e5593a84f9ac7" uuid = "ec8451be-7e33-11e9-00cf-bbf324bd1392" -version = "0.10.65" +version = "0.10.66" [[deps.LAME_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -849,9 +870,9 @@ version = "1.11.0" [[deps.LeftChildRightSiblingTrees]] deps = ["AbstractTrees"] -git-tree-sha1 = "fb6803dafae4a5d62ea5cab204b1e657d9737e7f" +git-tree-sha1 = "95ba48564903b43b2462318aa243ee79d81135ff" uuid = "1d6d02ad-be62-4b6b-8a6d-2f90e265016e" -version = "0.2.0" +version = "0.2.1" [[deps.LibCURL]] deps = ["LibCURL_jll", "MozillaCACerts_jll"] @@ -896,15 +917,15 @@ version = "1.18.0+0" [[deps.Libmount_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "a31572773ac1b745e0343fe5e2c8ddda7a37e997" +git-tree-sha1 = "3acf07f130a76f87c041cfb2ff7d7284ca67b072" uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9" -version = "2.41.0+0" +version = "2.41.2+0" [[deps.Libuuid_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "321ccef73a96ba828cd51f2ab5b9f917fa73945a" +git-tree-sha1 = "2a7a12fc0a4e7fb773450d17975322aa77142106" uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700" -version = "2.41.0+0" +version = "2.41.2+0" [[deps.LineSearches]] deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"] @@ -917,11 +938,23 @@ deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" version = "1.11.0" +[[deps.LinearMaps]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "7f6be2e4cdaaf558623d93113d6ddade7b916209" +uuid = "7a12625a-238d-50fd-b39a-03d52299707e" +version = "3.11.4" +weakdeps = ["ChainRulesCore", "SparseArrays", "Statistics"] + + [deps.LinearMaps.extensions] + LinearMapsChainRulesCoreExt = "ChainRulesCore" + LinearMapsSparseArraysExt = "SparseArrays" + LinearMapsStatisticsExt = "Statistics" + [[deps.LogDensityProblems]] deps = ["ArgCheck", "DocStringExtensions", "Random"] -git-tree-sha1 = "4e0128c1590d23a50dcdb106c7e2dbca99df85c0" +git-tree-sha1 = "d9625f27ded4ad726ceca7819394a4cc77ed25b3" uuid = "6fdf6af0-433a-55f7-b3ed-c6c6e0b8df7c" -version = "2.1.2" +version = "2.2.0" [[deps.LogExpFunctions]] deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] @@ -945,21 +978,43 @@ version = "1.11.0" [[deps.LoggingExtras]] deps = ["Dates", "Logging"] -git-tree-sha1 = "f02b56007b064fbfddb4c9cd60161b6dd0f40df3" +git-tree-sha1 = "f00544d95982ea270145636c181ceda21c4e2575" uuid = "e6f89c97-d47a-5376-807f-9c37f3926c36" -version = "1.1.0" +version = "1.2.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 = "e5afce7eaf5b5ca0d444bcb4dc4fd78c54cbbac0" +git-tree-sha1 = "a9fc7883eb9b5f04f46efb9a540833d1fad974b3" uuid = "bdcacae8-1622-11e9-2a5c-532679323890" -version = "0.12.172" -weakdeps = ["ChainRulesCore", "ForwardDiff", "SpecialFunctions"] +version = "0.12.173" [deps.LoopVectorization.extensions] ForwardDiffExt = ["ChainRulesCore", "ForwardDiff"] + ForwardDiffNNlibExt = ["ForwardDiff", "NNlib"] SpecialFunctionsExt = "SpecialFunctions" + [deps.LoopVectorization.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" + NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd" + SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" + +[[deps.LowRankApprox]] +deps = ["FFTW", "LinearAlgebra", "LowRankMatrices", "Nullables", "Random", "SparseArrays"] +git-tree-sha1 = "031af63ba945e23424815014ba0e59c28f5aed32" +uuid = "898213cb-b102-5a47-900c-97e73b919f73" +version = "0.5.5" + +[[deps.LowRankMatrices]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "59c5bb0708be6796604caec16d4357013dc3d132" +uuid = "e65ccdef-c354-471a-8090-89bec1c20ec3" +version = "1.0.2" +weakdeps = ["FillArrays"] + + [deps.LowRankMatrices.extensions] + LowRankMatricesFillArraysExt = "FillArrays" + [[deps.MCMCChains]] deps = ["AbstractMCMC", "AxisArrays", "Dates", "Distributions", "Formatting", "IteratorInterfaceExtensions", "KernelDensity", "LinearAlgebra", "MCMCDiagnosticTools", "MLJModelInterface", "NaturalSort", "OrderedCollections", "PrettyTables", "Random", "RecipesBase", "Serialization", "Statistics", "StatsBase", "StatsFuns", "TableTraits", "Tables"] git-tree-sha1 = "c659f7508035a7bdd5102aef2de028ab035f289a" @@ -974,15 +1029,15 @@ version = "0.2.1" [[deps.MKL_jll]] deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "oneTBB_jll"] -git-tree-sha1 = "5de60bc6cb3899cd318d80d627560fae2e2d99ae" +git-tree-sha1 = "282cadc186e7b2ae0eeadbd7a4dffed4196ae2aa" uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7" -version = "2025.0.1+1" +version = "2025.2.0+0" [[deps.MLJModelInterface]] -deps = ["REPL", "Random", "ScientificTypesBase", "StatisticalTraits"] -git-tree-sha1 = "66626f80d5807921045d539b4f7153b1d47c5f8a" +deps = ["InteractiveUtils", "REPL", "Random", "ScientificTypesBase", "StatisticalTraits"] +git-tree-sha1 = "ccaa3f7938890ee8042cc970ba275115428bd592" uuid = "e80e1ace-859a-464e-9ed9-23947d8ae3ea" -version = "1.11.1" +version = "1.12.0" [[deps.MacroTools]] git-tree-sha1 = "1e0228a030642014fe5cfe68c2c0a818f9e3f522" @@ -1007,9 +1062,9 @@ version = "0.1.2" [[deps.MathOptInterface]] deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "DataStructures", "ForwardDiff", "JSON3", "LinearAlgebra", "MutableArithmetics", "NaNMath", "OrderedCollections", "PrecompileTools", "Printf", "SparseArrays", "SpecialFunctions", "Test"] -git-tree-sha1 = "2f2c18c6acab9042557bdb0af8c3a14dd7b64413" +git-tree-sha1 = "a2cbab4256690aee457d136752c404e001f27768" uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" -version = "1.41.0" +version = "1.46.0" [[deps.MbedTLS_jll]] deps = ["Artifacts", "Libdl"] @@ -1038,9 +1093,9 @@ version = "2023.12.12" [[deps.MutableArithmetics]] deps = ["LinearAlgebra", "SparseArrays", "Test"] -git-tree-sha1 = "491bdcdc943fcbc4c005900d7463c9f216aabf4c" +git-tree-sha1 = "22df8573f8e7c593ac205455ca088989d0a2c7a0" uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" -version = "1.6.4" +version = "1.6.7" [[deps.NLSolversBase]] deps = ["ADTypes", "DifferentiationInterface", "Distributed", "FiniteDiff", "ForwardDiff"] @@ -1063,6 +1118,11 @@ version = "1.0.0" uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" version = "1.2.0" +[[deps.Nullables]] +git-tree-sha1 = "8f87854cc8f3685a60689d8edecaa29d2251979b" +uuid = "4d1e1d77-625e-5b40-9113-a560ec7a8ecd" +version = "1.0.0" + [[deps.OffsetArrays]] git-tree-sha1 = "117432e406b5c023f665fa73dc26e79ec3630151" uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" @@ -1096,15 +1156,15 @@ version = "0.8.5+0" [[deps.OpenSSH_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "OpenSSL_jll", "Zlib_jll"] -git-tree-sha1 = "cb7acd5d10aff809b4d0191dfe1956c2edf35800" +git-tree-sha1 = "301412a644646fdc0ad67d0a87487466b491e53d" uuid = "9bd350c2-7e96-507f-8002-3f2e150b4e1b" -version = "10.0.1+0" +version = "10.2.1+0" [[deps.OpenSSL_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "87510f7292a2b21aeff97912b0898f9553cc5c2c" +git-tree-sha1 = "f19301ae653233bc88b1810ae908194f07f8db9d" uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" -version = "3.5.1+0" +version = "3.5.4+0" [[deps.OpenSpecFun_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"] @@ -1197,9 +1257,9 @@ version = "1.2.1" [[deps.Preferences]] deps = ["TOML"] -git-tree-sha1 = "9306f6085165d270f7e3db02af26a400d580f5c6" +git-tree-sha1 = "0f27480397253da18fe2c12a4ba4eb9eb208bf3d" uuid = "21216c6a-2e73-6563-6e65-726566657250" -version = "1.4.3" +version = "1.5.0" [[deps.PrettyTables]] deps = ["Crayons", "LaTeXStrings", "Markdown", "PrecompileTools", "Printf", "Reexport", "StringManipulation", "Tables"] @@ -1230,9 +1290,9 @@ version = "0.1.5" [[deps.ProgressMeter]] deps = ["Distributed", "Printf"] -git-tree-sha1 = "13c5103482a8ed1536a54c08d0e742ae3dca2d42" +git-tree-sha1 = "fbb92c6c56b34e1a2c4c36058f68f332bec840e7" uuid = "92933f4c-e287-5a05-a399-4b506db050ca" -version = "1.10.4" +version = "1.11.0" [[deps.PtrArrays]] git-tree-sha1 = "1d36ef11a9aaf1e8b74dacc6a731dd1de8fd493d" @@ -1332,10 +1392,10 @@ uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f" version = "0.4.3+0" [[deps.SCS]] -deps = ["MathOptInterface", "PrecompileTools", "SCS_jll", "SparseArrays"] -git-tree-sha1 = "4aed85dec0209d638c241c34160999eaaf07965a" +deps = ["LinearAlgebra", "MathOptInterface", "OpenBLAS32_jll", "PrecompileTools", "SCS_jll", "SparseArrays"] +git-tree-sha1 = "48ec3c39787bc7b278789b9af17c157ea8774dae" uuid = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13" -version = "2.1.0" +version = "2.4.0" [deps.SCS.extensions] SCSSCS_GPU_jllExt = ["SCS_GPU_jll"] @@ -1346,10 +1406,10 @@ version = "2.1.0" SCS_MKL_jll = "3f2553a9-4106-52be-b7dd-865123654657" [[deps.SCS_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl", "OpenBLAS32_jll"] -git-tree-sha1 = "902cc4e42ecca21bbd74babf899b2a5b12add323" +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl", "libblastrampoline_jll"] +git-tree-sha1 = "05d6e31efa3debae6618dabee35dbd53cf4539d8" uuid = "f4f2fc5b-1d94-523c-97ea-2ab488bedf4b" -version = "3.2.7+0" +version = "300.200.900+0" [[deps.SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" @@ -1366,6 +1426,11 @@ git-tree-sha1 = "456f610ca2fbd1c14f5fcf31c6bfadc55e7d66e0" uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa" version = "0.6.43" +[[deps.SciMLPublic]] +git-tree-sha1 = "ed647f161e8b3f2973f24979ec074e8d084f1bee" +uuid = "431bcebd-1456-4ced-9d72-93c2757fff0b" +version = "1.0.0" + [[deps.ScientificTypesBase]] git-tree-sha1 = "a8e18eb383b5ecf1b5e6fc237eb39255044fd92b" uuid = "30f210dd-8aff-4c5f-94ba-8e64358c1161" @@ -1410,9 +1475,9 @@ version = "1.11.0" [[deps.SortingAlgorithms]] deps = ["DataStructures"] -git-tree-sha1 = "66e0a8e672a0bdfca2c3f5937efb8538b9ddc085" +git-tree-sha1 = "64d974c2e6fdf07f8155b5b2ca2ffa9069b608d9" uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c" -version = "1.2.1" +version = "1.2.2" [[deps.SparseArrays]] deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] @@ -1421,9 +1486,9 @@ version = "1.11.0" [[deps.SpecialFunctions]] deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] -git-tree-sha1 = "41852b8679f78c8d8961eeadc8f62cef861a52e3" +git-tree-sha1 = "f2685b435df2613e25fc10ad8c26dddb8640f547" uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "2.5.1" +version = "2.6.1" weakdeps = ["ChainRulesCore"] [deps.SpecialFunctions.extensions] @@ -1437,15 +1502,15 @@ version = "0.1.15" [[deps.StableRNGs]] deps = ["Random"] -git-tree-sha1 = "95af145932c2ed859b63329952ce8d633719f091" +git-tree-sha1 = "4f96c596b8c8258cc7d3b19797854d368f243ddc" uuid = "860ef19b-820b-49d6-a774-d7a799459cd3" -version = "1.0.3" +version = "1.0.4" [[deps.Static]] -deps = ["CommonWorldInvalidations", "IfElse", "PrecompileTools"] -git-tree-sha1 = "f737d444cb0ad07e61b3c1bef8eb91203c321eff" +deps = ["CommonWorldInvalidations", "IfElse", "PrecompileTools", "SciMLPublic"] +git-tree-sha1 = "49440414711eddc7227724ae6e570c7d5559a086" uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" -version = "1.2.0" +version = "1.3.1" [[deps.StaticArrayInterface]] deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "Static"] @@ -1460,9 +1525,9 @@ weakdeps = ["OffsetArrays", "StaticArrays"] [[deps.StaticArrays]] deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"] -git-tree-sha1 = "0feb6b9031bd5c51f9072393eb5ab3efd31bf9e4" +git-tree-sha1 = "b8693004b385c842357406e3af647701fe783f98" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.9.13" +version = "1.9.15" weakdeps = ["ChainRulesCore", "Statistics"] [deps.StaticArrays.extensions] @@ -1470,15 +1535,15 @@ weakdeps = ["ChainRulesCore", "Statistics"] StaticArraysStatisticsExt = "Statistics" [[deps.StaticArraysCore]] -git-tree-sha1 = "192954ef1208c7019899fbf8049e717f92959682" +git-tree-sha1 = "6ab403037779dae8c514bad259f32a447262455a" uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" -version = "1.4.3" +version = "1.4.4" [[deps.StatisticalTraits]] deps = ["ScientificTypesBase"] -git-tree-sha1 = "542d979f6e756f13f862aa00b224f04f9e445f11" +git-tree-sha1 = "89f86d9376acd18a1a4fbef66a56335a3a7633b8" uuid = "64bff920-2084-43da-a3e6-9bb72801c0c9" -version = "3.4.0" +version = "3.5.0" [[deps.Statistics]] deps = ["LinearAlgebra"] @@ -1520,6 +1585,20 @@ git-tree-sha1 = "159331b30e94d7b11379037feeb9b690950cace8" uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4" version = "1.11.0" +[[deps.StructUtils]] +deps = ["Dates", "UUIDs"] +git-tree-sha1 = "79529b493a44927dd5b13dde1c7ce957c2d049e4" +uuid = "ec057cc2-7a8d-4b58-b3b3-92acb9f63b42" +version = "2.6.0" + + [deps.StructUtils.extensions] + StructUtilsMeasurementsExt = ["Measurements"] + StructUtilsTablesExt = ["Tables"] + + [deps.StructUtils.weakdeps] + Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" + Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" + [[deps.StyledStrings]] uuid = "f489334b-da3d-4c2e-b8f0-e476e12c162b" version = "1.11.0" @@ -1643,9 +1722,9 @@ version = "1.11.0" [[deps.VectorizationBase]] deps = ["ArrayInterface", "CPUSummary", "HostCPUFeatures", "IfElse", "LayoutPointers", "Libdl", "LinearAlgebra", "SIMDTypes", "Static", "StaticArrayInterface"] -git-tree-sha1 = "4ab62a49f1d8d9548a1c8d1a75e5f55cf196f64e" +git-tree-sha1 = "d1d9a935a26c475ebffd54e9c7ad11627c43ea85" uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f" -version = "0.21.71" +version = "0.21.72" [[deps.VersionParsing]] git-tree-sha1 = "58d6e80b4ee071f5efd07fda82cb9fbe17200868" @@ -1713,15 +1792,15 @@ version = "0.2.7" [[deps.libaom_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "522c1df09d05a71785765d19c9524661234738e9" +git-tree-sha1 = "371cc681c00a3ccc3fbc5c0fb91f58ba9bec1ecf" uuid = "a4ae2306-e953-59d6-aa16-d00cac43593b" -version = "3.11.0+0" +version = "3.13.1+0" [[deps.libass_jll]] deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Zlib_jll"] -git-tree-sha1 = "e17c115d55c5fbb7e52ebedb427a0dca79d4484e" +git-tree-sha1 = "125eedcb0a4a0bba65b657251ce1d27c8714e9d6" uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0" -version = "0.15.2+0" +version = "0.17.4+0" [[deps.libblastrampoline_jll]] deps = ["Artifacts", "Libdl"] @@ -1752,10 +1831,10 @@ uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" version = "1.59.0+0" [[deps.oneTBB_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "d5a767a3bb77135a99e433afe0eb14cd7f6914c3" +deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl"] +git-tree-sha1 = "1350188a69a6e46f799d3945beef36435ed7262f" uuid = "1317d2d5-d96f-522e-a858-c73665f53c3e" -version = "2022.0.0+0" +version = "2022.0.0+1" [[deps.p7zip_jll]] deps = ["Artifacts", "Libdl"] @@ -1763,13 +1842,13 @@ uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" version = "17.4.0+2" [[deps.x264_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "4fea590b89e6ec504593146bf8b988b2c00922b2" +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "14cc7083fc6dff3cc44f2bc435ee96d06ed79aa7" uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a" -version = "2021.5.5+0" +version = "10164.0.1+0" [[deps.x265_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "ee567a171cce03570d77ad3a43e90218e38937a9" +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "e7b67590c14d487e734dcb925924c5dc43ec85f3" uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76" -version = "3.5.0+0" +version = "4.1.0+0" diff --git a/src/Emulator.jl b/src/Emulator.jl index b796ef508..a4caa6517 100644 --- a/src/Emulator.jl +++ b/src/Emulator.jl @@ -193,11 +193,7 @@ $(TYPEDSIGNATURES) Encode a new structure matrix in the input space (`"in"`) or output space (`"out"`). with the stored and initialized encoder schedule. """ -function encode_structure_matrix( - emulator::Emulator, - structure_mat::USorM, - in_or_out::AS, -) where {AS <: AbstractString, USorM <: Union{UniformScaling, AbstractMatrix}} +function encode_structure_matrix(emulator::Emulator, structure_mat, in_or_out::AS) where {AS <: AbstractString} return encode_with_schedule(get_encoder_schedule(emulator), structure_mat, in_or_out) end @@ -224,11 +220,7 @@ $(TYPEDSIGNATURES) Decode a new structure matrix in the input space (`"in"`) or output space (`"out"`). with the stored and initialized encoder schedule. """ -function decode_structure_matrix( - emulator::Emulator, - structure_mat::USorM, - in_or_out::AS, -) where {AS <: AbstractString, USorM <: Union{UniformScaling, AbstractMatrix}} +function decode_structure_matrix(emulator::Emulator, structure_mat, in_or_out::AS) where {AS <: AbstractString} return decode_with_schedule(get_encoder_schedule(emulator), structure_mat, in_or_out) end @@ -281,11 +273,11 @@ function predict( decoded_covariances = zeros(eltype(encoded_outputs), output_dim, output_dim, size(encoded_uncertainties)[end]) if var_or_cov == "var" for (i, col) in enumerate(eachcol(encoded_uncertainties)) - decoded_covariances[:, :, i] .= decode_structure_matrix(emulator, Diagonal(col), "out") + decoded_covariances[:, :, i] .= Matrix(decode_structure_matrix(emulator, Diagonal(col), "out")) end else # == "cov" for (i, mat) in enumerate(eachslice(encoded_uncertainties, dims = 3)) - decoded_covariances[:, :, i] .= decode_structure_matrix(emulator, mat, "out") + decoded_covariances[:, :, i] .= Matrix(decode_structure_matrix(emulator, mat, "out")) end end diff --git a/src/GaussianProcess.jl b/src/GaussianProcess.jl index 7f267216f..4f8a089b7 100644 --- a/src/GaussianProcess.jl +++ b/src/GaussianProcess.jl @@ -194,13 +194,7 @@ function build_models!( regularization = if isempty(output_structure_mats) 1.0 * ones(N_models) else - output_structure_mat = get_structure_mat(output_structure_mats) - if isa(output_structure_mat, UniformScaling) - output_structure_mat.λ * ones(N_models) - else - diag(output_structure_mat) - end - + output_structure_mat = diag(Matrix(get_structure_mat(output_structure_mats))) end regularization_noise = regularization .* gp.alg_reg_noise logstd_regularization_noise = log.(sqrt.(regularization_noise)) @@ -329,7 +323,7 @@ function build_models!( regularization = if isempty(output_structure_mats) 1.0 * ones(N_models) else - output_structure_mat = get_structure_mat(output_structure_mats) + output_structure_mat = Matrix(get_structure_mat(output_structure_mats)) if isa(output_structure_mat, UniformScaling) output_structure_mat.λ * ones(N_models) else @@ -435,7 +429,7 @@ AbstractGP currently does not (yet) learn hyperparameters internally. The follow regularization = if isempty(output_structure_mats) 1.0 * ones(N_models) else - output_structure_mat = get_structure_mat(output_structure_mats) + output_structure_mat = Matrix(get_structure_mat(output_structure_mats)) if isa(output_structure_mat, UniformScaling) output_structure_mat.λ * ones(N_models) else diff --git a/src/MarkovChainMonteCarlo.jl b/src/MarkovChainMonteCarlo.jl index 31d374040..06319d33a 100644 --- a/src/MarkovChainMonteCarlo.jl +++ b/src/MarkovChainMonteCarlo.jl @@ -33,6 +33,7 @@ export EmulatorPosteriorModel, accept_ratio, optimize_stepsize, get_posterior, + get_sample_kwargs, sample, esjd @@ -518,6 +519,13 @@ struct MCMCWrapper{AVV <: AbstractVector, AV <: AbstractVector} sample_kwargs::NamedTuple end +""" +$(TYPEDSIGNATURES) + +gets the NameTuple of keywords that are passed into the Sampler algorithm +""" +get_sample_kwargs(mcmc::MCMCWrapper) = mcmc.sample_kwargs + """ $(DocStringExtensions.TYPEDSIGNATURES) @@ -549,11 +557,15 @@ function MCMCWrapper( observation::AMorAV, prior::ParameterDistribution, em::Emulator; - init_params::AV, + init_params::AV = [], burnin::Int = 0, kwargs..., ) where {AV <: AbstractVector, AMorAV <: Union{AbstractVector, AbstractMatrix}} + if length(init_params) == 0 + init_params = mean(prior) + end + # make into iterable over vectors obs_slice = if observation isa AbstractVector{<:AbstractVector} observation @@ -587,6 +599,26 @@ function MCMCWrapper( return MCMCWrapper(prior, obs_slice, encoded_obs, log_posterior_map, mh_proposal_sampler, sample_kwargs) end +function MCMCWrapper( + mcmc_alg::MCMCProtocol, + observation::OB, + prior::PD, + em::EM; + kwargs..., +) where {OB <: Observation, PD <: ParameterDistribution, EM <: Emulator} + return MCMCWrapper(mcmc_alg, get_obs(observation), prior, em; kwargs...) +end + +function MCMCWrapper( + mcmc_alg::MCMCProtocol, + observation_series::OS, + prior::PD, + em::EM; + kwargs..., +) where {OS <: ObservationSeries, PD <: ParameterDistribution, EM <: Emulator} + observations = [get_obs(ob) for ob in get_observations(observation_series)] + return MCMCWrapper(mcmc_alg, observations, prior, em; kwargs...) +end """ $(DocStringExtensions.FUNCTIONNAME)([rng,] mcmc::MCMCWrapper, args...; kwargs...) diff --git a/src/ScalarRandomFeature.jl b/src/ScalarRandomFeature.jl index 0599b09a0..a6264a83f 100644 --- a/src/ScalarRandomFeature.jl +++ b/src/ScalarRandomFeature.jl @@ -404,12 +404,7 @@ function build_models!( regularization = if isempty(output_structure_mats) 1.0 * I(n_rfms) else - output_structure_mat = get_structure_mat(output_structure_mats) - if isa(output_structure_mat, UniformScaling) - output_structure_mat - else - Diagonal(output_structure_mat) - end + output_structure_mat = Diagonal(Matrix(get_structure_mat(output_structure_mats))) end @info ( diff --git a/src/Utilities.jl b/src/Utilities.jl index 298926c93..95a4b3824 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -1,18 +1,23 @@ module Utilities - - using DocStringExtensions using LinearAlgebra +using LinearMaps using Statistics using StatsBase using Random using ..EnsembleKalmanProcesses -EnsembleKalmanProcess = EnsembleKalmanProcesses.EnsembleKalmanProcess +const EKP = EnsembleKalmanProcesses +using EnsembleKalmanProcesses.ParameterDistributions using ..DataContainers -export get_training_points +using LowRankApprox +using TSVD +import LinearAlgebra: norm + +export get_training_points +export create_compact_linear_map export PairedDataContainerProcessor, DataContainerProcessor export create_encoder_schedule, initialize_and_encode_with_schedule!, @@ -21,10 +26,14 @@ export create_encoder_schedule, encode_data, encode_structure_matrix, decode_data, - decode_structure_matrix + decode_structure_matrix, + norm, + norm_linear_map, + isequal_linear, + encoder_kwargs_from -const StructureMatrix = Union{UniformScaling, AbstractMatrix} +const StructureMatrix = Union{UniformScaling, AbstractMatrix, AbstractVector, LinearMap} # The vector appears due to possible block-structured matrices (build=false) const StructureVector = Union{AbstractVector, AbstractMatrix} # In case of a matrix, the columns should be seen as vectors """ @@ -38,7 +47,7 @@ Extract the training points needed to train the Gaussian process regression. """ function get_training_points( - ekp::EnsembleKalmanProcess{FT, IT, P}, + ekp::EKP.EnsembleKalmanProcess{FT, IT, P}, train_iterations::Union{IT, AbstractVector{IT}}, ) where {FT, IT, P} @@ -63,20 +72,262 @@ function get_training_points( return training_points end +# Using Observation Objects: + +""" +$(TYPEDSIGNATURES) + +Extracts the relevant encoder kwargs from the observation as a NamedTuple. Contains, +- `:obs_noise_cov` as (unbuilt) noise covariance +- `:observation` as obs vector +""" +function encoder_kwargs_from(obs::OB) where {OB <: Observation} + return (; obs_noise_cov = get_obs_noise_cov(obs, build = false), observation = get_obs(obs)) +end + +""" +$(TYPEDSIGNATURES) + +Extracts the relevant encoder kwargs from the ObservationSeries as a NamedTuple. Assumes the same noise covariance for all observation vectors. Contains, +- `:obs_noise_cov` as (unbuilt) noise covariance of FIRST observation +- `:observation` as obs vector from all observations +""" +function encoder_kwargs_from(os::OS) where {OS <: ObservationSeries} + observations = get_observations(os) + obs_vec = [get_obs(obs) for obs in observations] + obs_noise_cov = get_obs_noise_cov(observations[1], build = false) + if !all([get_obs_noise_cov(observations[i], build = false) == obs_noise_cov for i in length(observations)]) + @warn(""" + Detected that observation covariances vary for different observations. + Encoder kwarg `:obs_noise_cov` will be set to the FIRST of these covariances for the purpose of data processing. + """) + end + return (; obs_noise_cov = obs_noise_cov, observation = obs_vec) +end + +""" +$(TYPEDSIGNATURES) + +Extracts the relevant encoder kwargs from the ParameterDistribution prior. Contains, +- `:prior_cov` as prior covariance +""" +function encoder_kwargs_from(prior::PD) where {PD <: ParameterDistribution} + return (; prior_cov = cov(prior)) +end + +## multiplication with observation covariance objects without building +""" +$(TYPEDSIGNATURES) + +Produces a linear map of type `LinearMap` that can evaluates the stacked actions of the structure matrix in compact form. by calling say `linear_map.f(x)` or `linear_map.fc(x)` to obtain `Ax`, or `A'x`. This particular type can be used by packages like `TSVD.jl` or `IterativeSolvers.jl` for further computations. + +This compact map constructs the following form of the Linear map f: + +1. get compact form svd-plus-d form "USVt + D" of the `blocks` +2. create the f via stacking `A.U * A.S * A.Vt * xblock + A.D * xblock for (A,xblock) in (As, x)` + +kwargs: +------ +When computing the svd internally from an abstract matrix +- `svd_dim_max=3000`: this switches to an approximate svd approach when applying to covariance matrices above dimension 3000 +- `psvd_or_tsvd="psvd"`: use psvd or tsvd for approximating svd for large matrices +- `tsvd_max_rank=50`: when using tsvd, what max rank to use. high rank = higher accuracy +- `psvd_kwargs=(; rtol=1e-2)`: when using psvd, what kwargs to pass. lower rtol = higher accuracy + +Recommended: quick & inaccurate -> slow and more accurate +- very large matrices - start with tsvd with very low rank, and increase +- mid-size matrices - psvd with very high rtol, and decrease +""" +function create_compact_linear_map( + A; + svd_dim_max = 3000, + psvd_or_tsvd = "psvd", + tsvd_max_rank = 50, + psvd_kwargs = (; rtol = 1e-1), +) + Avec = isa(A, AbstractVector) ? A : [A] + + # explicitly write the loop here: + Us = [] + Ss = [] + VTs = [] + ds = [] + batches = [] + shift = 0 + for a in Avec + bsize = 0 + if isa(a, UniformScaling) + throw( + ArgumentError( + "Detected `UniformScaling` (i.e. \"λI\") StructureMatrix, and unable to infer dimensionality. \n Please recast this as a diagonal matrix, defining \"λI(d)\" for dimension d", + ), + ) + end + if isa(a, AbstractMatrix) + if size(a, 1) <= svd_dim_max + svda = svd(a) + push!(Us, svda.U) + push!(Ss, svda.S) + push!(VTs, svda.Vt) + push!(ds, zeros(size(a, 1))) + bsize = size(a, 1) + else + if psvd_or_tsvd == "psvd" + svda = psvd(a; psvd_kwargs...) + push!(Us, svda[1]) + push!(Ss, svda[2]) + push!(VTs, svda[3]') + else + svda = tsvd_mat(a, min(tsvd_max_rank, size(a, 1) - 1)) + push!(Us, svda.U) + push!(Ss, svda.S) + push!(VTs, svda.Vt) + end + push!(ds, zeros(size(a, 1))) + bsize = size(a, 1) + end + elseif isa(a, SVD) + svda = a + push!(Us, svda.U) + push!(Ss, svda.S) + push!(VTs, svda.Vt) + push!(ds, zeros(size(a.U, 1))) + bsize = size(a.U, 1) + elseif isa(a, SVDplusD) + svda = a.svd_cov + diaga = (a.diag_cov).diag + push!(Us, svda.U) + push!(Ss, svda.S) + push!(VTs, svda.Vt) + push!(ds, diaga) + bsize = length(diaga) + end + + batch = (shift + 1):(shift + bsize) + push!(batches, batch) + shift = batch[end] + end + + # then create the LinearMap with entries (f(x) = A*x, f(x) = A'*x, size(A,1), size(A,2)) + # LinearMaps can only be applied to vectors in general, so we only provide this argumentation + + Amap = LinearMap( + x -> reduce( + vcat, + [U * (S .* (Vt * x[batch])) + d .* x[batch] for (U, S, Vt, d, batch) in zip(Us, Ss, VTs, ds, batches)], + ), + x -> reduce( + vcat, + [Vt' * (S .* (U' * x[batch])) + d .* x[batch] for (U, S, Vt, d, batch) in zip(Us, Ss, VTs, ds, batches)], + ), + sum(size(U, 1) for U in Us), + sum(size(Vt, 2) for Vt in VTs), + ) + + return Amap + +end + +""" +$(TYPEDSIGNATURES) + +Approximately computes the norm of a `LinearMap` object. For `Amap` associated with matrix `A`, `norm_linear_map(Amap,p)≈norm(A,p)`. Can be aliased as `norm()` + +kwargs +------ +- n_eval(=nothing): number of mat-vec products to apply in the approximation (larger is more accurate). default performs `size(map,2)` products +- rng(=Random.default_rng()): random number generator + +""" +function norm_linear_map(A::LM, p::Real = 2; n_eval = nothing, rng = Random.default_rng()) where {LM <: LinearMap} + m, n = size(A) + + # use unit-normalized gaussian vectors + n_basis = isa(n_eval, Nothing) ? n : n_eval + samples = randn(n, n_basis) + for i in 1:size(samples, 2) + samples[:, i] /= norm(samples[:, i]) + end + out = zeros(m, n_basis) + mul!(out, A, samples)# must use mul! for multiply with lin map to return a matrix) + norm_val = (n / n_basis)^(1 / p) * norm(out, p) + + return norm_val +end + +LinearAlgebra.norm(A::LM, p::Real = 2; lm_kwargs...) where {LM <: LinearMap} = norm_linear_map(A, p; lm_kwargs...) # Data processing tooling: abstract type DataProcessor end abstract type PairedDataContainerProcessor <: DataProcessor end # tools that operate on inputs and outputs abstract type DataContainerProcessor <: DataProcessor end # tools that operate on only inputs or outputs - # define how to have equality -Base.:(==)(a::DCP, b::DCP) where {DCP <: DataContainerProcessor} = - all(getfield(a, f) == getfield(b, f) for f in fieldnames(DCP)) -Base.:(==)(a::PDCP, b::PDCP) where {PDCP <: PairedDataContainerProcessor} = - all(getfield(a, f) == getfield(b, f) for f in fieldnames(PDCP)) +# this gets messy with LinearMaps, + +""" +$(TYPEDSIGNATURES) + +Tests equality for a LinearMap on a standard basis of the input space. Note that this operation requires a matrix multiply per input dimension so can be expensive. + +Kwargs: +------- +- n_eval (=nothing): the number of basis vectors to compare against (randomly selected without replacement if `n_eval < size(A,1)`) +- tol (=2*eps()): the tolerance for equality on evaluation per entry +- rng (=default_rng()): When provided, and `n_eval < size(A,1)`; a random subset of the basis is compared, using this `rng`. +- up_to_sign(=false): Only assess equality up to a sign-error (sufficient for e.g. encoder/decoder matrices) +""" +function isequal_linear( + A::LM1, + B::LM2; + tol = 2 * eps(), + n_eval = nothing, + rng = Random.default_rng(), + up_to_sign = false, +) where {LM1 <: LinearMap, LM2 <: LinearMap} + m, n = size(A) + if !(n == size(B, 2)) + @warn "Comparing equality of linear maps with size ($(m), $(n)) and ($(size(B,1)), $(size(B,2))). Was this intended?" + return false + end + if !(m == size(B, 1)) + @warn "Comparing equality of linear maps with size ($(m), $(n)) and ($(size(B,1)),$(size(B,2))). Was this intended?" + return false + end + + # test on standard basis (up to n_eval tests) + basis_id = isa(n_eval, Nothing) ? collect(1:n) : randperm(rng, n)[1:n_eval] + + e = vec(zeros(eltype(A), n)) + for j in basis_id + e[j] += 1 + if up_to_sign + AmB = abs.(A * e) - abs.(B * e) + else + AmB = A * e - B * e + end + if !(norm(AmB) <= sqrt(n) * tol) + return false + end + e[j] -= 1 + end + return true +end +function Base.:(==)(a::LM1, b::LM2) where {LM1 <: LinearMap, LM2 <: LinearMap} + n = size(a, 2) + if n < 1e4 # gets expensive + return isequal_linear(a, b; tol = 1e-12) + else + return isequal_linear(a, b; n_eval = Int(floor(sqrt(n))), tol = 1e-12) # 1e4 compares ~ 100 evals, 1e7 compares ~ 3000 evals + end +end +Base.:(==)(a::DCP1, b::DCP2) where {DCP1 <: DataContainerProcessor, DCP2 <: DataContainerProcessor} = + all(getfield(a, f) == getfield(b, f) for f in fieldnames(DCP1)) + +Base.:(==)(a::PDCP1, b::PDCP2) where {PDCP1 <: PairedDataContainerProcessor, PDCP2 <: PairedDataContainerProcessor} = + all(getfield(a, f) == getfield(b, f) for f in fieldnames(PDCP1)) #### function get_structure_vec(structure_vecs, name = nothing) @@ -197,7 +448,6 @@ This function creates the encoder scheduler that is also machine readable. E.g., enc_schedule = [ (DataProcessor1(...), "in"), (DataProcessor2(...), "out"), - (DataProcessor2(...), "out"), (PairedDataProcessor3(...),"in"), (DataProcessor4(...), "in"), (DataProcessor4(...), "out"), @@ -250,15 +500,32 @@ function initialize_and_encode_with_schedule!( ) where {VV <: AbstractVector, PDC <: PairedDataContainer} processed_io_pairs = deepcopy(io_pairs) + # We additionally convert the `mats` into a linear-maps for flexible handling of massive covariances. In a matrix-free manner input_structure_mats = deepcopy(input_structure_mats) if !isnothing(prior_cov) (input_structure_mats[:prior_cov] = prior_cov) end + for (key, val) in input_structure_mats + if isa(val, UniformScaling) # remove this annoying case immediately + input_dim = size(get_inputs(io_pairs), 1) + input_structure_mats[key] = create_compact_linear_map(val(input_dim)) + else + input_structure_mats[key] = create_compact_linear_map(val) + end + end output_structure_mats = deepcopy(output_structure_mats) if !isnothing(obs_noise_cov) (output_structure_mats[:obs_noise_cov] = obs_noise_cov) end + for (key, val) in output_structure_mats + if isa(val, UniformScaling) # remove this annoying case immediately + output_dim = size(get_outputs(io_pairs), 1) + output_structure_mats[key] = create_compact_linear_map(val(output_dim)) + else + output_structure_mats[key] = create_compact_linear_map(val) + end + end input_structure_vecs = deepcopy(input_structure_vecs) if !isnothing(prior_samples_in) @@ -342,9 +609,9 @@ Takes in an already initialized encoder schedule, and encodes a structure matrix """ function encode_with_schedule( encoder_schedule::VV, - structure_matrix::USorM, + structure_matrix::SM, in_or_out::AS, -) where {VV <: AbstractVector, USorM <: Union{UniformScaling, AbstractMatrix}, AS <: AbstractString} +) where {VV <: AbstractVector, SM <: StructureMatrix, AS <: AbstractString} if in_or_out ∉ ["in", "out"] bad_in_or_out(in_or_out) end @@ -368,14 +635,9 @@ Takes in an already initialized encoder schedule, and decodes a `DataContainer`, function decode_with_schedule( encoder_schedule::VV, io_pairs::PDC, - input_structure_mat::USorM1, - output_structure_mat::USorM2, -) where { - VV <: AbstractVector, - USorM1 <: Union{UniformScaling, AbstractMatrix}, - USorM2 <: Union{UniformScaling, AbstractMatrix}, - PDC <: PairedDataContainer, -} + input_structure_mat::SM1, + output_structure_mat::SM2, +) where {VV <: AbstractVector, SM1 <: StructureMatrix, SM2 <: StructureMatrix, PDC <: PairedDataContainer} processed_io_pairs = deepcopy(io_pairs) processed_input_structure_mat = deepcopy(input_structure_mat) processed_output_structure_mat = deepcopy(output_structure_mat) @@ -433,9 +695,9 @@ Takes in an already initialized encoder schedule, and decodes a structure matrix """ function decode_with_schedule( encoder_schedule::VV, - structure_matrix::USorM, + structure_matrix::SM, in_or_out::AS, -) where {VV <: AbstractVector, USorM <: Union{UniformScaling, AbstractMatrix}, AS <: AbstractString} +) where {VV <: AbstractVector, SM <: StructureMatrix, AS <: AbstractString} if in_or_out ∉ ["in", "out"] bad_in_or_out(in_or_out) end diff --git a/src/Utilities/canonical_correlation.jl b/src/Utilities/canonical_correlation.jl index c7ab6a47c..9e52cc1fd 100644 --- a/src/Utilities/canonical_correlation.jl +++ b/src/Utilities/canonical_correlation.jl @@ -186,7 +186,9 @@ Apply the `CanonicalCorrelation` encoder, on a columns-are-data matrix function encode_data(cc::CanonicalCorrelation, data::MM) where {MM <: AbstractMatrix} data_mean = get_data_mean(cc)[1] encoder_mat = get_encoder_mat(cc)[1] - return encoder_mat * (data .- data_mean) + out = zeros(size(encoder_mat, 1), size(data, 2)) + mul!(out, encoder_mat, data .- data_mean) + return out end """ @@ -197,7 +199,9 @@ Apply the `CanonicalCorrelation` decoder, on a columns-are-data matrix function decode_data(cc::CanonicalCorrelation, data::MM) where {MM <: AbstractMatrix} data_mean = get_data_mean(cc)[1] decoder_mat = get_decoder_mat(cc)[1] - return decoder_mat * data .+ data_mean + out = zeros(size(decoder_mat, 1), size(data, 2)) + mul!(out, decoder_mat, data) + return out .+ data_mean end """ diff --git a/src/Utilities/decorrelator.jl b/src/Utilities/decorrelator.jl index 9f237a827..aa94266f2 100644 --- a/src/Utilities/decorrelator.jl +++ b/src/Utilities/decorrelator.jl @@ -1,7 +1,14 @@ # included in Utilities.jl export Decorrelator, decorrelate_sample_cov, decorrelate_structure_mat, decorrelate -export get_data_mean, get_encoder_mat, get_decoder_mat, get_retain_var, get_decorrelate_with +export get_data_mean, + get_encoder_mat, + get_decoder_mat, + get_retain_var, + get_decorrelate_with, + get_n_totvar_samples, + get_max_rank, + get_psvd_kwargs """ $(TYPEDEF) @@ -22,10 +29,21 @@ The SVD is taken over the estimated covariance of the data. The data samples wil For `decorrelate(;decorrelate_with="combined")` (default): The SVD is taken to be the sum of structure matrix and estimated covariance. This may be more robust to ill-specification of structure matrix, or poor estimation of the sample covariance. +Depending on the size of the matrix, we perform different options of SVD: + +Small Matrix (dim < 3000): + use LinearAlgebra.svd(Matrix) +Large Matrix (dim > 3000): + if retain_var = 1.0 + use LowRankApprox.psvd(LinearMap; psvd_kwargs...) + if retain_var < 1.0 + use TSVD.tsvd(LinearMap) + + # Fields $(TYPEDFIELDS) """ -struct Decorrelator{VV1, VV2, VV3, FT, AS <: AbstractString} <: DataContainerProcessor +struct Decorrelator{VV1, VV2, VV3, FT, NT <: NamedTuple, AS <: AbstractString} <: DataContainerProcessor "storage for the data mean" data_mean::VV1 "the matrix used to perform encoding" @@ -34,6 +52,12 @@ struct Decorrelator{VV1, VV2, VV3, FT, AS <: AbstractString} <: DataContainerPro decoder_mat::VV3 "the fraction of variance to be retained after truncating singular values (1 implies no truncation)" retain_var::FT + "when retain_var < 1, number of samples to estimate the total variance. Larger values reduce the error in approximation at the cost of additional matrix-vector products." + n_totvar_samples::Int + "maximum dimension of subspace for `retain_var < 1`. The search may become expensive at large ranks, and therefore can be cut-off in this way" + max_rank::Int + "when `retain_var = 1`, the `psvd` algorithm from `LowRankApprox.jl` is used to decorrelate the space. here, kwargs can be passed in as a NamedTuple" + psvd_kwargs::NT "Switch to choose what form of matrix to use to decorrelate the data" decorrelate_with::AS "When given, use the structure matrix by this name if `decorrelate_with` uses structure matrices. When `nothing`, try to use the only present structure matrix instead." @@ -50,8 +74,24 @@ Constructs the `Decorrelator` struct. Users can add optional keyword arguments: - `"sample_cov"`, see [`decorrelate_sample_cov`](@ref) - `"combined"`, sums the `"sample_cov"` and `"structure_mat"` matrices """ -decorrelate(; retain_var::FT = Float64(1.0), decorrelate_with = "combined", structure_mat_name = nothing) where {FT} = - Decorrelator([], [], [], clamp(retain_var, FT(0), FT(1)), decorrelate_with, structure_mat_name) +decorrelate(; + retain_var::FT = Float64(1.0), + decorrelate_with = "combined", + structure_mat_name = nothing, + n_totvar_samples::Int = 500, + max_rank::Int = 100, + psvd_kwargs = (; rtol = 1e-5), +) where {FT} = Decorrelator( + [], + [], + [], + clamp(retain_var, FT(0), FT(1)), + n_totvar_samples, + max_rank, + psvd_kwargs, + decorrelate_with, + structure_mat_name, +) """ $(TYPEDSIGNATURES) @@ -59,8 +99,22 @@ $(TYPEDSIGNATURES) Constructs the `Decorrelator` struct, setting decorrelate_with = "sample_cov". Encoding data with this will ensure that the distribution of data samples after encoding will be `Normal(0,I)`. One can additionally add keywords: - `retain_var`[=`1.0`]: to project onto the leading singular vectors such that `retain_var` variance is retained """ -decorrelate_sample_cov(; retain_var::FT = Float64(1.0)) where {FT} = - Decorrelator([], [], [], clamp(retain_var, FT(0), FT(1)), "sample_cov", nothing) +decorrelate_sample_cov(; + retain_var::FT = Float64(1.0), + n_totvar_samples::Int = 500, + max_rank::Int = 100, + psvd_kwargs = (; rtol = 1e-5), +) where {FT} = Decorrelator( + [], + [], + [], + clamp(retain_var, FT(0), FT(1)), + n_totvar_samples, + max_rank, + psvd_kwargs, + "sample_cov", + nothing, +) """ $(TYPEDSIGNATURES) @@ -68,8 +122,23 @@ $(TYPEDSIGNATURES) Constructs the `Decorrelator` struct, setting decorrelate_with = "structure_mat". This encoding will transform a provided structure matrix into `I`. One can additionally add keywords: - `retain_var`[=`1.0`]: to project onto the leading singular vectors such that `retain_var` variance is retained """ -decorrelate_structure_mat(; retain_var::FT = Float64(1.0), structure_mat_name = nothing) where {FT} = - Decorrelator([], [], [], clamp(retain_var, FT(0), FT(1)), "structure_mat", structure_mat_name) +decorrelate_structure_mat(; + retain_var::FT = Float64(1.0), + structure_mat_name = nothing, + n_totvar_samples::Int = 500, + max_rank::Int = 100, + psvd_kwargs = (; rtol = 1e-3), +) where {FT} = Decorrelator( + [], + [], + [], + clamp(retain_var, FT(0), FT(1)), + n_totvar_samples, + max_rank, + psvd_kwargs, + "structure_mat", + structure_mat_name, +) """ $(TYPEDSIGNATURES) @@ -102,6 +171,27 @@ get_retain_var(dd::Decorrelator) = dd.retain_var """ $(TYPEDSIGNATURES) +returns the `n_totvar_samples` field of the `Decorrelator`. +""" +get_n_totvar_samples(dd::Decorrelator) = dd.n_totvar_samples + +""" +$(TYPEDSIGNATURES) + +returns the `max_rank` field of the `Decorrelator`. +""" +get_max_rank(dd::Decorrelator) = dd.max_rank + +""" +$(TYPEDSIGNATURES) + +returns the `psvd_kwargs` field of the `Decorrelator`. +""" +get_psvd_kwargs(dd::Decorrelator) = dd.psvd_kwargs + +""" +$(TYPEDSIGNATURES) + returns the `decorrelate_with` field of the `Decorrelator`. """ get_decorrelate_with(dd::Decorrelator) = dd.decorrelate_with @@ -123,36 +213,30 @@ Computes and populates the `data_mean` and `encoder_mat` and `decoder_mat` field function initialize_processor!( dd::Decorrelator, data::MM, - structure_matrices::Dict{Symbol, <:StructureMatrix}, - ::Dict{Symbol, <:StructureVector}, -) where {MM <: AbstractMatrix} + structure_matrices::Dict{Symbol, SM}, + ::Dict{Symbol, SV}, +) where {MM <: AbstractMatrix, SM <: StructureMatrix, SV <: StructureVector} if length(get_data_mean(dd)) == 0 push!(get_data_mean(dd), vec(mean(data, dims = 2))) end if length(get_encoder_mat(dd)) == 0 - # Can do tsvd here for large matrices + # To unify the actions here, we use LinearMaps.jl + # Note: can only multiply Linear Maps with vectors, not matrices... decorrelate_with = get_decorrelate_with(dd) if decorrelate_with == "structure_mat" - structure_matrix = get_structure_mat(structure_matrices, dd.structure_mat_name) - if isa(structure_matrix, UniformScaling) - data_dim = size(data, 1) - svdA = svd(structure_matrix(data_dim)) - rk = data_dim - else - svdA = svd(structure_matrix) - rk = rank(structure_matrix) - end + decorrelation_map = get_structure_mat(structure_matrices, dd.structure_mat_name) + elseif decorrelate_with == "sample_cov" - cd = cov(data, dims = 2) - svdA = svd(cd) - rk = rank(cd) + decorrelation_action = tsvd_cov_from_samples(data) + decorrelation_map = create_compact_linear_map(decorrelation_action) + elseif decorrelate_with == "combined" - structure_matrix = get_structure_mat(structure_matrices, dd.structure_mat_name) - spluscd = structure_matrix + cov(data, dims = 2) - svdA = svd(spluscd) - rk = rank(spluscd) + map1 = get_structure_mat(structure_matrices, dd.structure_mat_name) + decorrelation_action2 = tsvd_cov_from_samples(data) + map2 = create_compact_linear_map(decorrelation_action2) + decorrelation_map = map1 + map2 # x -> dm.maps[1].f(x) + dm.maps[2].f(x) else throw( ArgumentError( @@ -160,26 +244,108 @@ function initialize_processor!( ), ) end + + n_totvar_samples = get_n_totvar_samples(dd) + max_rank = get_max_rank(dd) + psvd_kwargs = get_psvd_kwargs(dd) ret_var = get_retain_var(dd) - if ret_var < 1.0 - sv_cumsum = cumsum(svdA.S) / sum(svdA.S) # variance contributions are (sing_val) for these matrices - trunc_val = minimum(findall(x -> (x > ret_var), sv_cumsum)) - @info " truncating at $(trunc_val)/$(length(sv_cumsum)) retaining $(100.0*sv_cumsum[trunc_val])% of the variance of the structure matrix" + + max_svd_size = 3000 # max size to perform true SVD + n_norm_compute = 10 # number of repeated norm computations to estimate total variance (if ret_var < 1) + S = [] + Vt = [] + + # Three methods of performing: svd, partial svd or truncated svd. + # svd: convert to matrix and perform actual SVD. [good for small matrix (d < ~3000)] + # psvd: approximates SVD to a given error tolerance (can truncate too) [good for full svd approx] + # tsvd: truncates to dimension where resulting SVD has the same s.vals as the true [good of low rank approx] + + if maximum(size(decorrelation_map)) < max_svd_size + # svd options + svdA = svd(Matrix(decorrelation_map)) + rk = rank(Matrix(decorrelation_map)) + + ret_var = get_retain_var(dd) + if ret_var < 1.0 + sv_cumsum = cumsum(svdA.S) / sum(svdA.S) # variance contributions are (sing_val) for these matrices + trunc_val = minimum(findall(x -> (x > ret_var), sv_cumsum)) + @info " truncating at $(trunc_val)/$(length(sv_cumsum)) retaining $(100.0*sv_cumsum[trunc_val])% of the variance of the structure matrix" + else + trunc_val = rk + if rk < size(data, 1) + @info " truncating at $(trunc_val)/$(size(data,1)), as low-rank data detected" + end + end + + push!(S, svdA.S[1:trunc_val]) + push!(Vt, svdA.Vt[1:trunc_val, :]) + else - trunc_val = rk - if rk < size(data, 1) - @info " truncating at $(trunc_val)/$(size(data,1)), as low-rank data detected" + if ret_var < 1.0 + # tsvd option + + norm_sq_approx = zeros(n_norm_compute) + for i in 1:n_norm_compute + # approximate frobenius norm^2 (= sum of s.v.^2 = total variance) + norm_sq_approx[i] = norm_linear_map(decorrelation_map, n_eval = n_totvar_samples)^2 + end + @info "relative error of total variance $(std(norm_sq_approx)/mean(norm_sq_approx))" + retain_var_max = 1 - std(norm_sq_approx) / mean(norm_sq_approx) + + ret_var_tmp = min(ret_var, retain_var_max) + + for rk in 1:max_rank + _, Stmp, Vtmp = tsvd(decorrelation_map, rk) + retain_var_current = sum(Stmp .^ 2) / (retain_var_max * mean(norm_sq_approx)) + if retain_var_current > ret_var_tmp + push!(S, Stmp) + push!(Vt, Vtmp') + @info " truncating at $(rk)/$(size(data,1)) retaining $(retain_var_current*100)% (+/-$(100*std(norm_sq_approx)/mean(norm_sq_approx)))% of the variance of the structure matrix" + break + end + if rk == max_rank + @warn "Maximum tolerance of SVD truncation hit ($(rank_max)), consider increasing `rank_max` in decorrelator, this may also result from `retain_var` being close to 1, while `n_totvar_samples` is low. Proceeding with final iterant." + push!(S, Stmp) + push!(Vt, Vtmp') + end + end + else + # psvd option + _, Stmp, Vtmp = psvd(decorrelation_map; psvd_kwargs...) # randomized algorithm to avoid building the matrix. + + if length(Stmp) < size(data, 1) + @info " truncating at $(length(Stmp))/$(size(data,1)), as low-rank data detected" + end + push!(S, Stmp) + push!(Vt, Vtmp') + end end - sqrt_inv_sv = Diagonal(1.0 ./ sqrt.(svdA.S[1:trunc_val])) - sqrt_sv = Diagonal(sqrt.(svdA.S[1:trunc_val])) - # as we have svd of cov-matrix we can use U or Vt - encoder_mat = sqrt_inv_sv * svdA.Vt[1:trunc_val, :] - decoder_mat = svdA.Vt[1:trunc_val, :]' * sqrt_sv + # Enforce a sign convention for the singular vectors (rows Vt have positive first entry) + for (j, row) in enumerate(eachrow(Vt[1])) + if row[1] < 0 + Vt[1][j, :] *= -1 + end + end - push!(get_encoder_mat(dd), encoder_mat) - push!(get_decoder_mat(dd), decoder_mat) + # we explicitly make the encoder/decoder maps + encoder_map = LinearMap( + x -> Diagonal(1.0 ./ sqrt.(S[1])) * Vt[1] * x, # Ax + x -> Vt[1]' * Diagonal(1.0 ./ sqrt.(S[1])) * x, # A'x + size(Vt[1], 1), # size(A,1) + size(Vt[1], 2), # size(A,2) + ) + + decoder_map = LinearMap( + x -> Vt[1]' * Diagonal(sqrt.(S[1])) * x, # Ax + x -> Diagonal(sqrt.(S[1])) * Vt[1] * x, # A'x + size(Vt[1]', 1), # size(A,1) + size(Vt[1]', 2), # size(A,2) + ) + + push!(get_encoder_mat(dd), encoder_map) + push!(get_decoder_mat(dd), decoder_map) end end @@ -192,7 +358,9 @@ Apply the `Decorrelator` encoder, on a columns-are-data matrix function encode_data(dd::Decorrelator, data::MM) where {MM <: AbstractMatrix} data_mean = get_data_mean(dd)[1] encoder_mat = get_encoder_mat(dd)[1] - return encoder_mat * (data .- data_mean) + out = zeros(size(encoder_mat, 1), size(data, 2)) + mul!(out, encoder_mat, data .- data_mean) + return out end """ @@ -203,13 +371,15 @@ Apply the `Decorrelator` decoder, on a columns-are-data matrix function decode_data(dd::Decorrelator, data::MM) where {MM <: AbstractMatrix} data_mean = get_data_mean(dd)[1] decoder_mat = get_decoder_mat(dd)[1] - return decoder_mat * data .+ data_mean + out = zeros(size(decoder_mat, 1), size(data, 2)) + mul!(out, decoder_mat, data) + return out .+ data_mean end """ $(TYPEDSIGNATURES) -Apply the `Decorrelator` encoder to a provided structure matrix +Apply the `Decorrelator` encoder to a provided structure matrix. If the structure matrix is a LinearMap, then the encoded structure matrix remains a LinearMap. """ function encode_structure_matrix(dd::Decorrelator, structure_matrix::SM) where {SM <: StructureMatrix} encoder_mat = get_encoder_mat(dd)[1] @@ -219,7 +389,7 @@ end """ $(TYPEDSIGNATURES) -Apply the `Decorrelator` decoder to a provided structure matrix +Apply the `Decorrelator` decoder to a provided structure matrix. If the structure matrix is a LinearMap, then the encoded structure matrix remains a LinearMap. """ function decode_structure_matrix(dd::Decorrelator, enc_structure_matrix::SM) where {SM <: StructureMatrix} decoder_mat = get_decoder_mat(dd)[1] diff --git a/src/Utilities/elementwise_scaler.jl b/src/Utilities/elementwise_scaler.jl index 6cca8a9d2..a6521cd8d 100644 --- a/src/Utilities/elementwise_scaler.jl +++ b/src/Utilities/elementwise_scaler.jl @@ -2,7 +2,8 @@ export UnivariateAffineScaling, ElementwiseScaler, QuartileScaling, MinMaxScaling, ZScoreScaling export quartile_scale, minmax_scale, zscore_scale -export get_type, get_shift, get_scale +export get_type, + get_shift, get_scale, get_data_encoder_mat, get_data_decoder_mat, get_struct_encoder_mat, get_struct_decoder_mat """ $(TYPEDEF) @@ -16,11 +17,26 @@ Different methods `T` will build different transformations: and are accessed with [`get_type`](@ref) """ -struct ElementwiseScaler{T, VV <: AbstractVector} <: DataContainerProcessor +struct ElementwiseScaler{ + T, + VV <: AbstractVector, + VV2 <: AbstractVector, + VV3 <: AbstractVector, + VV4 <: AbstractVector, + VV5 <: AbstractVector, +} <: DataContainerProcessor "storage for the shift applied to data" shift::VV "storage for the scaling" scale::VV + "the matrix used to perform encoding of data samples" + data_encoder_mat::VV2 + "the inverse of the the matrix used to perform encoding of data samples" + data_decoder_mat::VV3 + "the matrix used to perform encoding of structure_matrices" + struct_encoder_mat::VV4 + "the inverse of the the matrix used to perform encoding of structure_matrices" + struct_decoder_mat::VV5 end abstract type UnivariateAffineScaling end @@ -55,7 +71,7 @@ For multivariate standardization, see [`Decorrelator`](@ref) zscore_scale() = ElementwiseScaler(ZScoreScaling) ElementwiseScaler(::Type{UAS}) where {UAS <: UnivariateAffineScaling} = - ElementwiseScaler{UAS, Vector{Float64}}(Float64[], Float64[]) + ElementwiseScaler{UAS, Vector{Float64}, Vector, Vector, Vector, Vector}(Float64[], Float64[], [], [], [], []) """ $(TYPEDSIGNATURES) @@ -78,11 +94,40 @@ Gets the `scale` field of the `ElementwiseScaler` """ get_scale(es::ElementwiseScaler) = es.scale +""" +$(TYPEDSIGNATURES) + +returns the `data_encoder_mat` field of the `ElementwiseScaler`. +""" +get_data_encoder_mat(es::ElementwiseScaler) = es.data_encoder_mat + +""" +$(TYPEDSIGNATURES) + +returns the `data_decoder_mat` field of the `ElementwiseScaler`. +""" +get_data_decoder_mat(es::ElementwiseScaler) = es.data_decoder_mat + +""" +$(TYPEDSIGNATURES) + +returns the `struct_encoder_mat` field of the `ElementwiseScaler`. +""" +get_struct_encoder_mat(es::ElementwiseScaler) = es.struct_encoder_mat + +""" +$(TYPEDSIGNATURES) + +returns the `struct_decoder_mat` field of the `ElementwiseScaler`. +""" +get_struct_decoder_mat(es::ElementwiseScaler) = es.struct_decoder_mat + function Base.show(io::IO, es::ElementwiseScaler) out = "ElementwiseScaler: $(get_type(es))" print(io, out) end + function initialize_processor!( es::ElementwiseScaler, data::MM, @@ -90,6 +135,8 @@ function initialize_processor!( ) where {MM <: AbstractMatrix, QS <: QuartileScaling} quartiles_vec = [quantile(dd, [0.25, 0.5, 0.75]) for dd in eachrow(data)] quartiles_mat = reduce(hcat, quartiles_vec) # 3 rows: Q1, Q2, and Q3 + + # we use these more for saving readable data append!(get_shift(es), quartiles_mat[2, :]) append!(get_scale(es), (quartiles_mat[3, :] - quartiles_mat[1, :])) end @@ -120,6 +167,38 @@ function initialize_processor!(es::ElementwiseScaler, data::MM) where {MM <: Abs if length(get_shift(es)) == 0 T = get_type(es) initialize_processor!(es, data, T) + + # we explicitly make the encoder/decoder maps + data_encoder_map = LinearMap( + x -> (x .- get_shift(es)) ./ get_scale(es), # Ax + x -> (x .- get_shift(es)) ./ get_scale(es), # A'x + size(data, 1), # size(A,1) + size(data, 1), # size(A,2) + ) + data_decoder_map = LinearMap( + x -> x .* get_scale(es) .+ get_shift(es), # Ax + x -> x .* get_scale(es) .+ get_shift(es), # A'x + size(data, 1), # size(A,1) + size(data, 1), # size(A,2) + ) + # the encoder for the structure matrix does not have a shift + struct_encoder_map = LinearMap( + x -> x ./ get_scale(es), # Ax + x -> x ./ get_scale(es), # A'x + size(data, 1), # size(A,1) + size(data, 1), # size(A,2) + ) + struct_decoder_map = LinearMap( + x -> x .* get_scale(es), # Ax + x -> x .* get_scale(es), # A'x + size(data, 1), # size(A,1) + size(data, 1), # size(A,2) + ) + push!(get_data_encoder_mat(es), data_encoder_map) + push!(get_data_decoder_mat(es), data_decoder_map) + push!(get_struct_encoder_mat(es), struct_encoder_map) + push!(get_struct_decoder_mat(es), struct_decoder_map) + end end @@ -129,11 +208,9 @@ $(TYPEDSIGNATURES) Apply the `ElementwiseScaler` encoder, on a columns-are-data matrix """ function encode_data(es::ElementwiseScaler, data::MM) where {MM <: AbstractMatrix} - out = deepcopy(data) - for i in 1:size(out, 1) - out[i, :] .-= get_shift(es)[i] - out[i, :] /= get_scale(es)[i] - end + out = zeros(size(data)) + enc = get_data_encoder_mat(es)[1] + mul!(out, enc, data) # must use this form to get matrix output of enc*out return out end @@ -143,11 +220,9 @@ $(TYPEDSIGNATURES) Apply the `ElementwiseScaler` decoder, on a columns-are-data matrix """ function decode_data(es::ElementwiseScaler, data::MM) where {MM <: AbstractMatrix} - out = deepcopy(data) - for i in 1:size(out, 1) - out[i, :] *= get_scale(es)[i] - out[i, :] .+= get_shift(es)[i] - end + out = zeros(size(data)) + dec = get_data_decoder_mat(es)[1] + mul!(out, dec, data) # must use this form to get matrix output of dec*out return out end @@ -167,17 +242,19 @@ initialize_processor!( """ $(TYPEDSIGNATURES) -Apply the `ElementwiseScaler` encoder to a provided structure matrix +Apply the `ElementwiseScaler` encoder to a provided structure matrix. If the structure matrix is a LinearMap, then the encoded structure matrix remains a LinearMap. """ function encode_structure_matrix(es::ElementwiseScaler, structure_matrix::SM) where {SM <: StructureMatrix} - return Diagonal(1 ./ get_scale(es)) * structure_matrix * Diagonal(1 ./ get_scale(es)) + encoder_mat = get_struct_encoder_mat(es)[1] + return encoder_mat * structure_matrix * encoder_mat' end """ $(TYPEDSIGNATURES) -Apply the `ElementwiseScaler` decoder to a provided structure matrix +Apply the `ElementwiseScaler` decoder to a provided structure matrix. If the structure matrix is a LinearMap, then the encoded structure matrix remains a LinearMap. """ function decode_structure_matrix(es::ElementwiseScaler, enc_structure_matrix::SM) where {SM <: StructureMatrix} - return Diagonal(get_scale(es)) * enc_structure_matrix * Diagonal(get_scale(es)) + decoder_mat = get_struct_decoder_mat(es)[1] + return decoder_mat * enc_structure_matrix * decoder_mat' end diff --git a/src/VectorRandomFeature.jl b/src/VectorRandomFeature.jl index 6d5d4f5a8..05eeccb2f 100644 --- a/src/VectorRandomFeature.jl +++ b/src/VectorRandomFeature.jl @@ -461,7 +461,7 @@ function build_models!( regularization = if isempty(output_structure_mats) 1.0 * I else - output_structure_mat = get_structure_mat(output_structure_mats) + output_structure_mat = Matrix(get_structure_mat(output_structure_mats)) if !isposdef(output_structure_mat) println("RF output structure matrix is not positive definite, correcting for use as a regularizer") posdef_correct(output_structure_mat) diff --git a/test/Emulator/runtests.jl b/test/Emulator/runtests.jl index 5eec9a6c0..898972c97 100644 --- a/test/Emulator/runtests.jl +++ b/test/Emulator/runtests.jl @@ -10,7 +10,7 @@ using LinearAlgebra using CalibrateEmulateSample.Emulators using CalibrateEmulateSample.DataContainers using CalibrateEmulateSample.Utilities - +using CalibrateEmulateSample.EnsembleKalmanProcesses #build an unknown type struct MLTester <: Emulators.MachineLearningTool end @@ -45,12 +45,11 @@ struct MLTester <: Emulators.MachineLearningTool end default_encoder = (decorrelate_sample_cov(), "in_and_out") # for these inputs this is the default enc_sch = create_encoder_schedule(default_encoder) enc_io_pairs, enc_I_in, enc_I_out = initialize_and_encode_with_schedule!(enc_sch, io_pairs; obs_noise_cov = 1.0 * I) - @test get_encoder_schedule(em)[1][1] == enc_sch[1][1] # inputs: proc - @test get_encoder_schedule(em)[1][2] == enc_sch[1][2] # inputs: apply_to - @test get_encoder_schedule(em)[2][1] == enc_sch[2][1] # outputs... - @test get_encoder_schedule(em)[2][2] == enc_sch[2][2] - @test get_data(get_encoded_io_pairs(em)) == get_data(enc_io_pairs) - @test get_data(get_encoded_io_pairs(em)) == get_data(enc_io_pairs) + # NB this gives encoder up to sign + tol = 1e-12 + @test get_encoder_schedule(em) == enc_sch # inputs: proc + @test all(isapprox.(get_inputs(get_encoded_io_pairs(em)), get_inputs(enc_io_pairs), atol = tol)) + @test all(isapprox.(get_outputs(get_encoded_io_pairs(em)), get_outputs(enc_io_pairs), atol = tol)) @test isempty(enc_I_in) #NB - encoders all tested in Utilities here just testing some API @@ -70,8 +69,8 @@ struct MLTester <: Emulators.MachineLearningTool end @test isapprox(norm(decoded_I - 1.0 * I), 0, atol = tol * d * d) # test obs_noise_cov (check the warning at the start) - @test_logs (:warn,) (:info,) (:info,) (:warn,) Emulator(gp, io_pairs, obs_noise_cov = Σ) - @test_logs (:warn,) (:info,) (:info,) (:warn,) Emulator( + @test_logs (:warn,) (:info,) (:warn,) (:info,) (:warn,) Emulator(gp, io_pairs, obs_noise_cov = Σ) + @test_logs (:warn,) (:info,) (:warn,) (:info,) (:warn,) Emulator( gp, io_pairs; encoder_kwargs = (; prior_cov = 4.0 * I, obs_noise_cov = 2.0 * I), diff --git a/test/MarkovChainMonteCarlo/runtests.jl b/test/MarkovChainMonteCarlo/runtests.jl index 7546c8718..3775939eb 100644 --- a/test/MarkovChainMonteCarlo/runtests.jl +++ b/test/MarkovChainMonteCarlo/runtests.jl @@ -4,6 +4,7 @@ using Distributions using GaussianProcesses using Test +using CalibrateEmulateSample.EnsembleKalmanProcesses using CalibrateEmulateSample.MarkovChainMonteCarlo const MCMC = MarkovChainMonteCarlo using CalibrateEmulateSample.ParameterDistributions @@ -251,11 +252,15 @@ function mcmc_test_template( target_acc = 0.25, return_samples = false, ) - if !isa(obs_sample, AbstractVecOrMat) + if isa(obs_sample, Real) obs_sample = reshape(collect(obs_sample), 1) # scalar -> Vector end init_params = vec(collect(init_params)) # scalar or Vector -> Vector + mcmc = MCMCWrapper(mcmc_alg, obs_sample, prior, em) # without ICs + + @test all(isapprox.(getfield(get_sample_kwargs(mcmc), :initial_params), mean(prior))) + mcmc = MCMCWrapper(mcmc_alg, obs_sample, prior, em; init_params = init_params) # First let's run a short chain to determine a good step size new_step = optimize_stepsize(mcmc; init_stepsize = step, N = 5000, target_acc = target_acc) @@ -301,7 +306,7 @@ end # mcmc setup mcmc_params = Dict( :mcmc_alg => RWMHSampling(), - :obs_sample => obs_sample, + :obs_sample => Observation(Dict("samples" => obs_sample, "covariances" => 1.0 * I, "names" => "test")), :init_params => transform_constrained_to_unconstrained(prior, [2.0]), :step => 0.25, :rng => rng, @@ -415,30 +420,44 @@ end @info "Posterior mean: $(posterior_mean_2) ≈ $(mle)" # test with many slightly different samples - # as vec of vec - obs_sample2 = [obs_sample + 0.01 * randn(length(obs_sample)) for i in 1:100] - mcmc_params2 = deepcopy(mcmc_params) - mcmc_params2[:obs_sample] = obs_sample2 - mcmc_params2[:step] = 0.025 # less uncertainty -> smaller step - new_step, posterior_mean_1 = mcmc_test_template(prior, σ2_y, em_1; exp_name = "gpjl-samples", mcmc_params2...) - @test isapprox(new_step, 0.025; atol = 0.025) - # difference between mean_1 and ground truth comes from MCMC convergence and GP sampling - esjd2 = esjd(chain_2) - @info "ESJD (-many-samples)= $esjd2" - @info "Posterior mean: $(posterior_mean_1) ≈ $(mle)" - @test isapprox(posterior_mean_1, mle; atol = 2e-1) + obs_sample_vecs = [] + exp_vec_names = [] + # as a vec of vec + obs_sample2 = [obs_sample + 0.01 * randn(length(obs_sample)) for i in 1:100] + push!(exp_vec_names, "gpjl-samples-mat") + push!(obs_sample_vecs, obs_sample2) - # as column matrix + # as a column mat obs_sample2mat = reduce(hcat, obs_sample2) - mcmc_params2mat = deepcopy(mcmc_params) - mcmc_params2mat[:obs_sample] = obs_sample2mat - new_step, posterior_mean_1 = - mcmc_test_template(prior, σ2_y, em_1; exp_name = "gpjl-samples-mat", mcmc_params2mat...) - @test isapprox(new_step, 0.025; atol = 0.025) - # difference between mean_1 and ground truth comes from MCMC convergence and GP sampling - @info "Posterior mean: $(posterior_mean_1) ≈ $(mle)" - @test isapprox(posterior_mean_1, mle; atol = 2e-1) + push!(obs_sample_vecs, obs_sample2mat) + push!(exp_vec_names, "gpjl-samples") + + # as an observation series + obs_vec = [] + for i in 1:100 + ob = Observation( + Dict("samples" => obs_sample + 0.01 * randn(length(obs_sample)), "covariances" => I, "names" => "test"), + ) + push!(obs_vec, ob) + end + obs_series = ObservationSeries(obs_vec) + push!(obs_sample_vecs, obs_series) + push!(exp_vec_names, "gpjl-samples-series") + + # run tests... + for (os, exp_name) in zip(obs_sample_vecs, exp_vec_names) + mcmc_params2 = deepcopy(mcmc_params) + mcmc_params2[:obs_sample] = os + mcmc_params2[:step] = 0.025 # less uncertainty -> smaller step + new_step, posterior_mean_1 = mcmc_test_template(prior, σ2_y, em_1; exp_name = exp_name, mcmc_params2...) + @test isapprox(new_step, 0.025; atol = 0.025) + # difference between mean_1 and ground truth comes from MCMC convergence and GP sampling + esjd2 = esjd(chain_2) + @info "ESJD (vec:$(exp_name))= $esjd2" + @info "Posterior mean: $(posterior_mean_1) ≈ $(mle)" + @test isapprox(posterior_mean_1, mle; atol = 2e-1) + end # test with integer data obs_sample3 = [4] diff --git a/test/RandomFeature/runtests.jl b/test/RandomFeature/runtests.jl index 3778e8c25..52863fcaf 100644 --- a/test/RandomFeature/runtests.jl +++ b/test/RandomFeature/runtests.jl @@ -340,14 +340,14 @@ rng = Random.MersenneTwister(seed) @test_throws ArgumentError Emulator(srfi_bad, iopairs) # test under repeats - @test_logs (:info,) (:info,) (:warn,) Emulator( + @test_logs (:info,) (:warn,) (:info,) (:warn,) Emulator( srfi, iopairs; encoder_kwargs = (; obs_noise_cov = obs_noise_cov), ) Emulator(srfi, iopairs; encoder_kwargs = (; obs_noise_cov = obs_noise_cov)) @test length(get_rfms(srfi)) == n_srfi - @test_logs (:info,) (:info,) (:warn,) Emulator( + @test_logs (:info,) (:warn,) (:info,) (:warn,) Emulator( vrfi, iopairs; encoder_kwargs = (; obs_noise_cov = obs_noise_cov), diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index f235ed0f3..f9140fca0 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -7,6 +7,7 @@ using LinearAlgebra using CalibrateEmulateSample.Utilities using CalibrateEmulateSample.EnsembleKalmanProcesses using CalibrateEmulateSample.DataContainers +using CalibrateEmulateSample.ParameterDistributions @testset "Utilities" begin @@ -49,6 +50,94 @@ end # Seed for pseudo-random number generator rng = Random.MersenneTwister(4154) + # test getting encoder kwargs from observations and series + # from prior + pd = constrained_gaussian("name", 0.0, 2.0, -10, Inf, repeat = 5) + encoder_kwargs = encoder_kwargs_from(pd) + @test isapprox(norm(encoder_kwargs.prior_cov - cov(pd)), 0, atol = 1e-12) + + # from observation + osample = [1.0, 2.0] + ocov = 4.0 * I(2) + obs = Observation(Dict("samples" => osample, "covariances" => ocov, "names" => "test")) + encoder_kwargs = encoder_kwargs_from(obs) + @test all(isapprox.(norm.(encoder_kwargs.obs_noise_cov - get_obs_noise_cov(obs, build = false)), 0.0, atol = 1e-12)) + @test all(isapprox.(encoder_kwargs.observation - get_obs(obs), 0.0, atol = 1e-12)) + + # from observation series + obs2 = Observation(Dict("samples" => osample .+ 1.0, "covariances" => 2.0 .* ocov, "names" => "test")) + obs_series = ObservationSeries([obs, obs2]) + encoder_kwargs = encoder_kwargs_from(obs_series) + @test all(isapprox.(norm.(encoder_kwargs.obs_noise_cov - get_obs_noise_cov(obs, build = false)), 0, atol = 1e-12)) + diff = encoder_kwargs.observation - [get_obs(obs), get_obs(obs2)] + @test all([all(isapprox.(dd, 0.0, atol = 1e-12)) for dd in diff]) + + # Building a blocked noise covariance with different types to test creating linear maps + # Abstract mats (small and large (> svd_dim_max)), SVD, SVDplusD + d = 1000 + A = [] + svd_dim_max = 999 # defaul 4000, reduce here for speed/memory + X = randn(100, d) + X = X' * X # make posdef + push!(A, X[1:50, 1:50]) + push!(A, X) + + μ = zeros(d) + sam = rand(d, 30) + Σ = tsvd_cov_from_samples(sam) + push!(A, Σ) + + D = Diagonal(0.001 * sqrt.(collect(1:d))) + ΣpD = SVDplusD(Σ, D) + push!(A, ΣpD) + + + as = [50, d, d, d] + cas = cumsum(as) + block_id = [1:cas[1], (cas[1] + 1):cas[2], (cas[2] + 1):cas[3], (cas[3] + 1):cas[4]] + Afull = zeros(sum(as), sum(as)) + + Afull[1:cas[1], 1:cas[1]] = A[1] + Afull[(cas[1] + 1):cas[2], (cas[1] + 1):cas[2]] = A[2] + Afull[(cas[2] + 1):cas[3], (cas[2] + 1):cas[3]] = A[3].U * Diagonal(A[3].S) * A[3].Vt + Afull[(cas[3] + 1):cas[4], (cas[3] + 1):cas[4]] = + A[4].svd_cov.U * Diagonal(A[4].svd_cov.S) * A[4].svd_cov.Vt + A[4].diag_cov + + + @test_throws ArgumentError create_compact_linear_map(3 * I) # needs dims + + for svd_type in ["psvd", "tsvd"] + psvd_kwargs = (; rtol = 1e-3) # make very small for testing + tsvd_max_rank = 30 # often only stable small ranks + Amap = create_compact_linear_map( + A; + svd_dim_max = svd_dim_max, + psvd_or_tsvd = svd_type, + psvd_kwargs = psvd_kwargs, + tsvd_max_rank = tsvd_max_rank, + ) + Amat_from_map = Matrix(Amap) + + # some extra tests for the ==/is_equal overrides + @test Amap == Amap + Amap2 = create_compact_linear_map([A[1]]) + @test !(Amap == Amap2) + + for (i, ids) in enumerate(block_id) + if i == 2 # the case where we do psvd on a full matrix - can be very high error expected (particularly tsvd + if svd_type == "psvd" + @test norm(Afull[ids, ids] - Amat_from_map[ids, ids]) < 1e-14 * length(ids)^2 + else + @test norm(Afull[ids, ids] - Amat_from_map[ids, ids]) < 1e-2 * length(ids)^2 # bad! + end + + else + @test norm(Afull[ids, ids] - Amat_from_map[ids, ids]) < 1e-14 * length(ids)^2 + end + end + end + + # Tests for get_structure_vec and get_structure_mat structure_vecs = Dict("a" => [1, 2, 3], "b" => [4, 5, 6]) structure_mats = Dict("c" => [1 2; 3 4], "d" => [5 6; 7 8]) @@ -72,7 +161,7 @@ end zs = zscore_scale() mm = minmax_scale() qq = quartile_scale() - QQ = ElementwiseScaler{QuartileScaling, Vector{Int}}([1], [2]) + QQ = ElementwiseScaler{QuartileScaling, Vector{Int}, Vector, Vector, Vector, Vector}([1], [2], [3], [4], [5], [6]) @test isa(zs, ElementwiseScaler) @test get_type(zs) == ZScoreScaling @test isa(mm, ElementwiseScaler) @@ -81,7 +170,10 @@ end @test get_type(qq) == QuartileScaling @test get_shift(QQ) == [1] @test get_scale(QQ) == [2] - + @test get_data_encoder_mat(QQ) == [3] + @test get_data_decoder_mat(QQ) == [4] + @test get_struct_encoder_mat(QQ) == [5] + @test get_struct_decoder_mat(QQ) == [6] dd = decorrelate() @test get_retain_var(dd) == 1.0 @test get_decorrelate_with(dd) == "combined" @@ -91,10 +183,13 @@ end dd3 = decorrelate_structure_mat(retain_var = 0.7) @test get_retain_var(dd3) == 0.7 @test get_decorrelate_with(dd3) == "structure_mat" - DD = Decorrelator([1], [2], [3], 1.0, "test", nothing) + DD = Decorrelator([1], [2], [3], 1.0, 4, 5, (; test = 6), "test", nothing) @test get_data_mean(DD) == [1] @test get_encoder_mat(DD) == [2] @test get_decoder_mat(DD) == [3] + @test get_n_totvar_samples(DD) == 4 + @test get_max_rank(DD) == 5 + @test get_psvd_kwargs(DD) == (; test = 6) cc = canonical_correlation() @@ -154,7 +249,7 @@ end (canonical_correlation(retain_var = 0.95), "in_and_out"), ] - lossless = [fill(true, 6); fill(false, 4)] # are these lossy approximations? + lossless = [fill(true, 6); fill(false, 4)] # are these lossless approximations? # functional test pipeline tol = 1e-12 @@ -381,7 +476,82 @@ end @test_throws ArgumentError encode_with_schedule(encoder_schedule, prior_cov, "bad") @test_throws ArgumentError decode_with_schedule(encoder_schedule, encoded_pc, "bad") +end + + + + +@testset "Decorrelator: Large observational covariance" begin + + # loop over output dim + ds = [10, 10, 100, 1_000, 10_000, 100_000] # 1_000_000 + dts_struct = [] + dts_comb = [] + dts_samp = [] + for (d_idx, d) in enumerate(ds) + @info " " + @info "Testing decorrelating dimension: $d" + @info " " + #build some quick data + noise + m = 50 + p = 10 + x = rand(p, m) #R^3 + y = rand(d, m) #R^6 + + # "noise" + μ = zeros(d) + sam = rand(d, 30) + Σ = tsvd_cov_from_samples(sam) + D = Diagonal(0.001 * sqrt.(collect(1:d))) + noise_samples = Σ.U * (sqrt.(Σ.S) .* Σ.Vt) * rand(MvNormal(μ, I), m) + y += noise_samples + io_pairs = PairedDataContainer(x, y, data_are_columns = true) + + dt = @elapsed begin + enc1 = (decorrelate_sample_cov(), "in_and_out") # for these inputs this is the default + + enc_sch1 = create_encoder_schedule(enc1) + enc_io_pairs, enc_in, enc_out = + initialize_and_encode_with_schedule!(enc_sch1, io_pairs; obs_noise_cov = [Σ]) + end + push!(dts_samp, dt) + + dt = @elapsed begin + enc2 = [(decorrelate_sample_cov(), "in"), (decorrelate_structure_mat(retain_var = 0.95), "out")] # for these inputs this is the default + enc_sch2 = create_encoder_schedule(enc2) + enc_io_pairs, enc_in, enc_out = + initialize_and_encode_with_schedule!(enc_sch2, io_pairs; obs_noise_cov = [Σ]) + end + push!(dts_struct, dt) + + dt = @elapsed begin + enc3 = [(decorrelate_sample_cov(), "in"), (decorrelate(retain_var = 0.95), "out")] # for these inputs this is the default + enc_sch3 = create_encoder_schedule(enc3) + enc_io_pairs, enc_in, enc_out = + initialize_and_encode_with_schedule!(enc_sch3, io_pairs; obs_noise_cov = [Σ]) + end + push!(dts_comb, dt) + + end + + dts_tols = 2 * [dts_samp[2], dts_struct[2], dts_comb[2]] # baseline value + for i in 1:length(ds) + + if i == 1 + println("dimension decorr-sample decorr-structure decorr-combined") + else + + # increate the tolerances by the factor + tols = 5 * ds[i] / ds[2] * dts_tols # should be < 5 x scaling + println("$(ds[i]), $(dts_samp[i]), $(dts_struct[i]), $(dts_comb[i])") + if any([dts_samp[i], dts_struct[i], dts_comb[i]] .> tols) + @error("∟ timings have exceeded linear scaling") + end + + end + end + # @tests? end