From c7e83dc15343277d7a5a5a6d1efc273a9a56a82c Mon Sep 17 00:00:00 2001 From: Thomas Jackson Date: Sat, 4 Dec 2021 16:52:34 -0500 Subject: [PATCH] Commit .dev/climaformat.jl Auto-generate signatures in documentation Use DocStringExtensions to auto-generate signatures and typedefs, since manually maintained docs are at risk of going out of sync with the code. Auto-generate signatures in documentation Use DocStringExtensions to auto-generate signatures and typedefs, since manually maintained docs are at risk of going out of sync with the code. Update Manifest for build on Julia 1.6 Tweak docstring grammar/formatting Ensure docs for constructors are generated Run climaformat.jl --- .dev/Project.toml | 5 + .dev/clima_formatter_options.jl | 8 + .dev/climaformat.jl | 85 +++ .dev/hooks/pre-commit | 41 ++ Manifest.toml | 92 ++- docs/Manifest.toml | 100 ++-- docs/make.jl | 47 +- docs/src/API/MarkovChainMonteCarlo.md | 11 +- docs/src/API/Utilities.md | 6 + docs/src/installation_instructions.md | 17 +- .../Emulator/GaussianProcess/learn_noise.jl | 46 +- examples/Emulator/GaussianProcess/plot_GP.jl | 208 ++++--- examples/Lorenz/GModel.jl | 525 +++++++++--------- examples/Lorenz/Lorenz_example.jl | 255 +++++---- examples/ci/linkfig.jl | 4 +- examples/deprecated/Cloudy/Cloudy_example.jl | 105 ++-- examples/deprecated/Cloudy/GModel.jl | 77 ++- src/CalibrateEmulateSample.jl | 19 +- src/Emulator.jl | 241 ++++---- src/GaussianProcess.jl | 148 +++-- src/MCMC_sampler_defaults.jl | 251 +++++++++ src/MarkovChainMonteCarlo.jl | 148 +++-- src/Utilities.jl | 53 +- test/Emulator/runtests.jl | 135 ++--- test/GaussianProcess/runtests.jl | 140 +++-- test/MarkovChainMonteCarlo/runtests.jl | 204 ++++--- test/Utilities/runtests.jl | 38 +- test/runtests.jl | 12 +- 28 files changed, 1770 insertions(+), 1251 deletions(-) create mode 100644 .dev/Project.toml create mode 100644 .dev/clima_formatter_options.jl create mode 100644 .dev/climaformat.jl create mode 100644 .dev/hooks/pre-commit create mode 100644 src/MCMC_sampler_defaults.jl diff --git a/.dev/Project.toml b/.dev/Project.toml new file mode 100644 index 000000000..d86279dc7 --- /dev/null +++ b/.dev/Project.toml @@ -0,0 +1,5 @@ +[deps] +JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" + +[compat] +JuliaFormatter = "0.3" diff --git a/.dev/clima_formatter_options.jl b/.dev/clima_formatter_options.jl new file mode 100644 index 000000000..97926de3f --- /dev/null +++ b/.dev/clima_formatter_options.jl @@ -0,0 +1,8 @@ +clima_formatter_options = ( + indent = 4, + margin = 120, + always_for_in = true, + whitespace_typedefs = true, + whitespace_ops_in_indices = true, + remove_extra_newlines = false, +) diff --git a/.dev/climaformat.jl b/.dev/climaformat.jl new file mode 100644 index 000000000..713a5189d --- /dev/null +++ b/.dev/climaformat.jl @@ -0,0 +1,85 @@ +#!/usr/bin/env julia +# +# This is an adapted version of format.jl from JuliaFormatter with the +# following license: +# +# MIT License +# Copyright (c) 2019 Dominique Luna +# +# Permission is hereby granted, free of charge, to any person obtaining a +# copy of this software and associated documentation files (the +# "Software"), to deal in the Software without restriction, including +# without limitation the rights to use, copy, modify, merge, publish, +# distribute, sublicense, and/or sell copies of the Software, and to permit +# persons to whom the Software is furnished to do so, subject to the +# following conditions: +# +# The above copyright notice and this permission notice shall be included +# in all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN +# NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, +# DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR +# OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE +# USE OR OTHER DEALINGS IN THE SOFTWARE. +# +using Pkg +Pkg.activate(@__DIR__) +Pkg.instantiate() + +using JuliaFormatter + +include("clima_formatter_options.jl") + +help = """ +Usage: climaformat.jl [flags] [FILE/PATH]... + +Formats the given julia files using the CLIMA formatting options. If paths +are given it will format the julia files in the paths. Otherwise, it will +format all changed julia files. + + -v, --verbose + Print the name of the files being formatted with relevant details. + + -h, --help + Print this message. +""" + +function parse_opts!(args::Vector{String}) + i = 1 + opts = Dict{Symbol, Union{Int, Bool}}() + while i ≤ length(args) + arg = args[i] + if arg[1] != '-' + i += 1 + continue + end + if arg == "-v" || arg == "--verbose" + opt = :verbose + elseif arg == "-h" || arg == "--help" + opt = :help + else + error("invalid option $arg") + end + if opt in (:verbose, :help) + opts[opt] = true + deleteat!(args, i) + end + end + return opts +end + +opts = parse_opts!(ARGS) +if haskey(opts, :help) + write(stdout, help) + exit(0) +end +if isempty(ARGS) + filenames = readlines(`git ls-files "*.jl"`) +else + filenames = ARGS +end + +format(filenames; clima_formatter_options..., opts...) diff --git a/.dev/hooks/pre-commit b/.dev/hooks/pre-commit new file mode 100644 index 000000000..0e56ed080 --- /dev/null +++ b/.dev/hooks/pre-commit @@ -0,0 +1,41 @@ +#!/usr/bin/env julia +# +# Called by git-commit with no arguments. This checks to make sure that all +# .jl files are indented correctly before a commit is made. +# +# To enable this hook, make this file executable and copy it in +# $GIT_DIR/hooks. + +toplevel_directory = chomp(read(`git rev-parse --show-toplevel`, String)) + +using Pkg +Pkg.activate(joinpath(toplevel_directory, ".dev")) +Pkg.instantiate() + +using JuliaFormatter + +include(joinpath(toplevel_directory, ".dev", "clima_formatter_options.jl")) + +needs_format = false + +for diffoutput in split.(readlines(`git diff --name-status --cached`)) + status = diffoutput[1] + filename = diffoutput[end] + (!startswith(status, "D") && endswith(filename, ".jl")) || continue + + a = read(`git show :$filename`, String) + b = format_text(a; clima_formatter_options...) + + if a != b + fullfilename = joinpath(toplevel_directory, filename) + + @error """File $filename needs to be indented with: + julia $(joinpath(toplevel_directory, ".dev", "climaformat.jl")) $fullfilename + and added to the git index via + git add $fullfilename + """ + global needs_format = true + end +end + +exit(needs_format ? 1 : 0) diff --git a/Manifest.toml b/Manifest.toml index c34040746..3c20a1bab 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -22,9 +22,9 @@ uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" [[ArrayInterface]] deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"] -git-tree-sha1 = "6e8fada11bb015ecf9263f64b156f98b546918c7" +git-tree-sha1 = "c933ce606f6535a7c7b98e1d86d5d1014f730596" uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -version = "5.0.5" +version = "5.0.7" [[Artifacts]] uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" @@ -50,17 +50,11 @@ git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2" uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0" version = "1.0.8+0" -[[Calculus]] -deps = ["LinearAlgebra"] -git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad" -uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9" -version = "0.5.1" - [[ChainRulesCore]] deps = ["Compat", "LinearAlgebra", "SparseArrays"] -git-tree-sha1 = "c9a6160317d1abe9c44b3beb367fd448117679ca" +git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f" uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "1.13.0" +version = "1.14.0" [[ChangesOfVariables]] deps = ["ChainRulesCore", "LinearAlgebra", "Test"] @@ -88,9 +82,9 @@ version = "0.3.0" [[Compat]] deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"] -git-tree-sha1 = "96b0bc6c52df76506efc8a441c6cf1adcb1babc4" +git-tree-sha1 = "b153278a25dd42c65abbf4e62344f9d22e59191b" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "3.42.0" +version = "3.43.0" [[CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] @@ -114,15 +108,15 @@ uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f" version = "4.1.1" [[DataAPI]] -git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8" +git-tree-sha1 = "fb5f5316dd3fd4c5e7c30a24d50643b73e37cd40" uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" -version = "1.9.0" +version = "1.10.0" [[DataFrames]] deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"] -git-tree-sha1 = "ae02104e835f219b8930c7664b8012c93475c340" +git-tree-sha1 = "6c19003824cbebd804a51211fd3bbd81bf1ecad5" uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" -version = "1.3.2" +version = "1.3.3" [[DataStructures]] deps = ["Compat", "InteractiveUtils", "OrderedCollections"] @@ -151,9 +145,9 @@ version = "1.0.3" [[DiffRules]] deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"] -git-tree-sha1 = "dd933c4ef7b4c270aacd4eb88fa64c147492acf0" +git-tree-sha1 = "28d605d9a0ac17118fe2c5e9ce0fbb76c3ceb120" uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" -version = "1.10.0" +version = "1.11.0" [[Distances]] deps = ["LinearAlgebra", "SparseArrays", "Statistics", "StatsAPI"] @@ -181,12 +175,6 @@ version = "0.8.6" deps = ["ArgTools", "LibCURL", "NetworkOptions"] uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" -[[DualNumbers]] -deps = ["Calculus", "NaNMath", "SpecialFunctions"] -git-tree-sha1 = "90b158083179a6ccbce2c7eb1446d5bf9d7ae571" -uuid = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74" -version = "0.6.7" - [[ElasticArrays]] deps = ["Adapt"] git-tree-sha1 = "a0fcc1bb3c9ceaf07e1d0529c9806ce94be6adf9" @@ -201,9 +189,9 @@ version = "0.2.2" [[EnsembleKalmanProcesses]] deps = ["Convex", "Distributions", "DocStringExtensions", "LinearAlgebra", "Random", "SCS", "SparseArrays", "Statistics", "StatsBase"] -git-tree-sha1 = "e0c9d7714fff35dbeaa6bdc7bbaa2d03006bf073" +git-tree-sha1 = "5e7b26419a9693f8b9223b640f0d8bb08f7432b0" uuid = "aa8a2aa5-91d8-4396-bcef-d4f2ec43552d" -version = "0.5.2" +version = "0.5.3" [[FastGaussQuadrature]] deps = ["LinearAlgebra", "SpecialFunctions", "StaticArrays"] @@ -231,9 +219,9 @@ version = "0.4.2" [[ForwardDiff]] deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions", "StaticArrays"] -git-tree-sha1 = "1bd6fc0c344fc0cbee1f42f8d2e7ec8253dda2d2" +git-tree-sha1 = "34e6147e7686a101c245f12dba43b743c7afda96" uuid = "f6369f11-7733-5829-9624-2563aa707210" -version = "0.10.25" +version = "0.10.27" [[Future]] deps = ["Random"] @@ -245,12 +233,6 @@ git-tree-sha1 = "9cf8ba8037e332b1be14c71e549143e68c42a22d" uuid = "891a1506-143c-57d2-908e-e1f8e92e6de9" version = "0.12.4" -[[HypergeometricFunctions]] -deps = ["DualNumbers", "LinearAlgebra", "SpecialFunctions", "Test"] -git-tree-sha1 = "65e4589030ef3c44d3b90bdc5aac462b4bb05567" -uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a" -version = "0.3.8" - [[IfElse]] git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1" uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" @@ -300,9 +282,9 @@ version = "0.21.3" [[LDLFactorizations]] deps = ["AMD", "LinearAlgebra", "SparseArrays", "Test"] -git-tree-sha1 = "399bbe845e06e1c2d44ebb241f554d45eaf66788" +git-tree-sha1 = "736e01b9b2d443c4e3351aebe551b8a374ab9c05" uuid = "40e66cde-538c-5869-a4ad-c39174c6795b" -version = "0.8.1" +version = "0.8.2" [[LibCURL]] deps = ["LibCURL_jll", "MozillaCACerts_jll"] @@ -335,9 +317,9 @@ uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" [[LogExpFunctions]] deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"] -git-tree-sha1 = "58f25e56b706f95125dcb796f39e1fb01d913a71" +git-tree-sha1 = "44a7b7bb7dd1afe12bac119df6a7e540fa2c96bc" uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" -version = "0.3.10" +version = "0.3.13" [[Logging]] uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" @@ -433,9 +415,9 @@ version = "0.12.3" [[Parsers]] deps = ["Dates"] -git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7" +git-tree-sha1 = "1285416549ccfcdf0c50d4997a94331e88d68413" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "2.2.3" +version = "2.3.1" [[Pkg]] deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] @@ -443,9 +425,9 @@ uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" [[PooledArrays]] deps = ["DataAPI", "Future"] -git-tree-sha1 = "db3a23166af8aebf4db5ef87ac5b00d36eb771e2" +git-tree-sha1 = "28ef6c7ce353f0b35d0df0d5930e0d072c1f5b9b" uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720" -version = "1.4.0" +version = "1.4.1" [[PositiveFactorizations]] deps = ["LinearAlgebra"] @@ -455,9 +437,9 @@ version = "0.2.4" [[Preferences]] deps = ["TOML"] -git-tree-sha1 = "d3538e7f8a790dc8903519090857ef8e1283eecd" +git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d" uuid = "21216c6a-2e73-6563-6e65-726566657250" -version = "1.2.5" +version = "1.3.0" [[PrettyTables]] deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"] @@ -475,9 +457,9 @@ uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79" [[ProgressMeter]] deps = ["Distributed", "Printf"] -git-tree-sha1 = "afadeba63d90ff223a6a48d2009434ecee2ec9e8" +git-tree-sha1 = "d7a7aef8f8f2d537104f170139553b14dfe39fe9" uuid = "92933f4c-e287-5a05-a399-4b506db050ca" -version = "1.7.1" +version = "1.7.2" [[PyCall]] deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"] @@ -588,15 +570,15 @@ version = "1.8.4" [[Static]] deps = ["IfElse"] -git-tree-sha1 = "87e9954dfa33fd145694e42337bdd3d5b07021a6" +git-tree-sha1 = "b1f1f60bf4f25d8b374480fb78c7b9785edf95fd" uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" -version = "0.6.0" +version = "0.6.2" [[StaticArrays]] deps = ["LinearAlgebra", "Random", "Statistics"] -git-tree-sha1 = "74fb527333e72ada2dd9ef77d98e4991fb185f04" +git-tree-sha1 = "cd56bf18ed715e8b09f06ef8c6b781e6cdc49911" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.4.1" +version = "1.4.4" [[Statistics]] deps = ["LinearAlgebra", "SparseArrays"] @@ -604,9 +586,9 @@ uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" [[StatsAPI]] deps = ["LinearAlgebra"] -git-tree-sha1 = "c3d8ba7f3fa0625b062b82853a7d5229cb728b6b" +git-tree-sha1 = "c82aaa13b44ea00134f8c9c89819477bd3986ecd" uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" -version = "1.2.1" +version = "1.3.0" [[StatsBase]] deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"] @@ -615,10 +597,10 @@ uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" version = "0.33.16" [[StatsFuns]] -deps = ["ChainRulesCore", "HypergeometricFunctions", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"] -git-tree-sha1 = "25405d7016a47cf2bd6cd91e66f4de437fd54a07" +deps = ["ChainRulesCore", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"] +git-tree-sha1 = "5950925ff997ed6fb3e985dcce8eb1ba42a0bbe7" uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c" -version = "0.9.16" +version = "0.9.18" [[SuiteSparse]] deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] diff --git a/docs/Manifest.toml b/docs/Manifest.toml index d95d44d36..f08274459 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -1,99 +1,95 @@ # This file is machine-generated - editing it directly is not advised -[[Base64]] +julia_version = "1.7.1" +manifest_format = "2.0" + +[[deps.ANSIColoredPrinters]] +git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c" +uuid = "a4c015fc-c6ff-483c-b24f-f7ea428134e9" +version = "0.0.1" + +[[deps.Base64]] uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" -[[Dates]] +[[deps.Dates]] deps = ["Printf"] uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" -[[Distributed]] -deps = ["Random", "Serialization", "Sockets"] -uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" - -[[DocStringExtensions]] -deps = ["LibGit2", "Markdown", "Pkg", "Test"] -git-tree-sha1 = "9d4f64f79012636741cf01133158a54b24924c32" +[[deps.DocStringExtensions]] +deps = ["LibGit2"] +git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b" uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.8.4" +version = "0.8.6" -[[Documenter]] -deps = ["Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"] -git-tree-sha1 = "3ebb967819b284dc1e3c0422229b58a40a255649" +[[deps.Documenter]] +deps = ["ANSIColoredPrinters", "Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"] +git-tree-sha1 = "f425293f7e0acaf9144de6d731772de156676233" uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -version = "0.26.3" +version = "0.27.10" -[[IOCapture]] -deps = ["Logging"] -git-tree-sha1 = "377252859f740c217b936cebcd918a44f9b53b59" +[[deps.IOCapture]] +deps = ["Logging", "Random"] +git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a" uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89" -version = "0.1.1" +version = "0.2.2" -[[InteractiveUtils]] +[[deps.InteractiveUtils]] deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" -[[JSON]] +[[deps.JSON]] deps = ["Dates", "Mmap", "Parsers", "Unicode"] -git-tree-sha1 = "81690084b6198a2e1da36fcfda16eeca9f9f24e4" +git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37" uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" -version = "0.21.1" +version = "0.21.2" -[[LibGit2]] -deps = ["Printf"] +[[deps.LibGit2]] +deps = ["Base64", "NetworkOptions", "Printf", "SHA"] uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" -[[Libdl]] -uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" - -[[Logging]] +[[deps.Logging]] uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" -[[Markdown]] +[[deps.Markdown]] deps = ["Base64"] uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" -[[Mmap]] +[[deps.Mmap]] uuid = "a63ad114-7e13-5084-954f-fe012c677804" -[[Parsers]] +[[deps.NetworkOptions]] +uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" + +[[deps.Parsers]] deps = ["Dates"] -git-tree-sha1 = "c8abc88faa3f7a3950832ac5d6e690881590d6dc" +git-tree-sha1 = "d7fa6237da8004be601e19bd6666083056649918" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "1.1.0" +version = "2.1.3" -[[Pkg]] -deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"] -uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" - -[[Printf]] +[[deps.Printf]] deps = ["Unicode"] uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" -[[REPL]] -deps = ["InteractiveUtils", "Markdown", "Sockets"] +[[deps.REPL]] +deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" -[[Random]] -deps = ["Serialization"] +[[deps.Random]] +deps = ["SHA", "Serialization"] uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -[[SHA]] +[[deps.SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" -[[Serialization]] +[[deps.Serialization]] uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" -[[Sockets]] +[[deps.Sockets]] uuid = "6462fe0b-24de-5631-8697-dd941f90decc" -[[Test]] -deps = ["Distributed", "InteractiveUtils", "Logging", "Random"] +[[deps.Test]] +deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -[[UUIDs]] -deps = ["Random", "SHA"] -uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" - -[[Unicode]] +[[deps.Unicode]] uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" diff --git a/docs/make.jl b/docs/make.jl index c90af5b4c..b60b14770 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -6,14 +6,12 @@ using Documenter # using DocumenterCitations #---------- -api = ["CalibrateEmulateSample" => [ - "Emulators" => [ - "General Emulator" => "API/Emulators.md", - "Gaussian Process" => "API/GaussianProcess.md", - ], - "MarkovChainMonteCarlo" => "API/MarkovChainMonteCarlo.md", - "Utilities" => "API/Utilities.md" - ], +api = [ + "CalibrateEmulateSample" => [ + "Emulators" => ["General Emulator" => "API/Emulators.md", "Gaussian Process" => "API/GaussianProcess.md"], + "MarkovChainMonteCarlo" => "API/MarkovChainMonteCarlo.md", + "Utilities" => "API/Utilities.md", + ], ] examples = ["Lorenz example" => "examples/lorenz_example.md"] @@ -27,27 +25,24 @@ pages = [ #---------- -format = Documenter.HTML( - collapselevel = 1, - prettyurls = !isempty(get(ENV, "CI", "")) -) +format = Documenter.HTML(collapselevel = 1, prettyurls = !isempty(get(ENV, "CI", ""))) makedocs( - sitename = "CalibrateEmulateSample.jl", - authors = "CliMA Contributors", - format = format, - pages = pages, - modules = [CalibrateEmulateSample], - doctest = false, - strict = true, - clean = true, - checkdocs = :none, + sitename = "CalibrateEmulateSample.jl", + authors = "CliMA Contributors", + format = format, + pages = pages, + modules = [CalibrateEmulateSample], + doctest = false, + strict = true, + clean = true, + checkdocs = :none, ) if !isempty(get(ENV, "CI", "")) - deploydocs( - repo = "github.com/CliMA/CalibrateEmulateSample.jl.git", - versions = ["stable" => "v^", "v#.#.#", "dev" => "dev"], - push_preview = true, - ) + deploydocs( + repo = "github.com/CliMA/CalibrateEmulateSample.jl.git", + versions = ["stable" => "v^", "v#.#.#", "dev" => "dev"], + push_preview = true, + ) end diff --git a/docs/src/API/MarkovChainMonteCarlo.md b/docs/src/API/MarkovChainMonteCarlo.md index 846a94d7f..97174354b 100644 --- a/docs/src/API/MarkovChainMonteCarlo.md +++ b/docs/src/API/MarkovChainMonteCarlo.md @@ -1,9 +1,6 @@ -# Utilities +# MarkovChainMonteCarlo -```@meta -CurrentModule = CalibrateEmulateSample.MarkovChainMonteCarlo +```@autodocs +Modules = [CalibrateEmulateSample.MarkovChainMonteCarlo] +Order = [:module, :type, :function] ``` - -```@docs -MCMC -``` \ No newline at end of file diff --git a/docs/src/API/Utilities.md b/docs/src/API/Utilities.md index e69de29bb..71f2f3bab 100644 --- a/docs/src/API/Utilities.md +++ b/docs/src/API/Utilities.md @@ -0,0 +1,6 @@ +# Utilities + +```@autodocs +Modules = [CalibrateEmulateSample.Utilities] +Order = [:module, :type, :function] +``` diff --git a/docs/src/installation_instructions.md b/docs/src/installation_instructions.md index dc92f4496..421c118c7 100644 --- a/docs/src/installation_instructions.md +++ b/docs/src/installation_instructions.md @@ -2,7 +2,9 @@ ### Installing CalibrateEmulateSample.jl -Currently CalibrateEmulateSample (CES) depends on some external python dependencies including `scikit-learn` wrapped by ScikitLearn.jl which requires a couple extra installation steps: +Currently CalibrateEmulateSample (CES) depends on some external python dependencies +including `scikit-learn` wrapped by ScikitLearn.jl, which requires a couple extra +installation steps: First clone the project into a new local repository @@ -11,8 +13,8 @@ First clone the project into a new local repository > cd CalibrateEmulateSample.jl ``` -Install and build the project dependencies. -Given that CES depends on python packages it is easiest to set the project to use it's own +Install and build the project dependencies. Given that CES depends on python packages +it is easiest to set the project to use its own [Conda](https://docs.conda.io/en/latest/miniconda.html) environment variable (set by exporting the ENV variable `PYTHON=""`). @@ -21,10 +23,11 @@ Given that CES depends on python packages it is easiest to set the project to us ``` -The `scikit-learn` package then has to be installed if using a Julia project specific Conda environment: +The `scikit-learn` package then has to be installed if using a Julia project-specific +Conda environment: ``` -> PYTHON="" julia --project= -e 'using Conda; Conda.add("scikit-learn")' +> PYTHON="" julia --project -e 'using Conda; Conda.add("scikit-learn")' ``` @@ -43,7 +46,7 @@ You need to first build the top-level project before building the documentation: ``` cd CalibrateEmulateSample.jl -julia --project -e 'using Pkg; Pkg.instantiate() +julia --project -e 'using Pkg; Pkg.instantiate()' ``` Then you can build the project documentation under the `docs/` sub-project: @@ -53,4 +56,4 @@ julia --project=docs/ -e 'using Pkg; Pkg.instantiate()' julia --project=docs/ docs/make.jl ``` -The locally rendered HTML documentation can be viewed at `docs/build/index.html` \ No newline at end of file +The locally rendered HTML documentation can be viewed at `docs/build/index.html`. diff --git a/examples/Emulator/GaussianProcess/learn_noise.jl b/examples/Emulator/GaussianProcess/learn_noise.jl index d9a985161..f28c9583c 100644 --- a/examples/Emulator/GaussianProcess/learn_noise.jl +++ b/examples/Emulator/GaussianProcess/learn_noise.jl @@ -51,10 +51,10 @@ Random.seed!(rng_seed) gppackage = GPJL() prediction_type = YType() gauss_proc_1 = GaussianProcess( - gppackage; - kernel=nothing, # use default squared exponential kernel - prediction_type=prediction_type, - noise_learn=true + gppackage; + kernel = nothing, # use default squared exponential kernel + prediction_type = prediction_type, + noise_learn = true, ) # Generate training data (x-y pairs, where x ∈ ℝ ᵖ, y ∈ ℝ ᵈ) @@ -72,11 +72,11 @@ X = 2.0 * π * rand(p, n) # G(x1, x2) g1x = sin.(X[1, :]) .+ cos.(X[2, :]) g2x = sin.(X[1, :]) .- cos.(X[2, :]) -gx = zeros(2,n) -gx[1,:] = g1x -gx[2,:] = g2x +gx = zeros(2, n) +gx[1, :] = g1x +gx[2, :] = g2x # Add noise η -μ = zeros(d) +μ = zeros(d) Σ = 0.1 * [[0.8, 0.2] [0.2, 0.5]] # d x d noise_samples = rand(MvNormal(μ, Σ), n) @@ -93,16 +93,11 @@ Y = gx .+ noise_samples # in the training phase via an optimization procedure. # Because of the SVD transformation applied to the output, we expect the # learned noise to be close to 1. -iopairs = PairedDataContainer(X,Y,data_are_columns=true) +iopairs = PairedDataContainer(X, Y, data_are_columns = true) @assert get_inputs(iopairs) == X @assert get_outputs(iopairs) == Y -emulator_1 = Emulator( - gauss_proc_1, - iopairs, - obs_noise_cov=Σ, - normalize_inputs=true -) +emulator_1 = Emulator(gauss_proc_1, iopairs, obs_noise_cov = Σ, normalize_inputs = true) println("build GP with ", n, " training points") #optimize the hyperparameters to best fit the GP. optimize_hyperparameters!(emulator_1) @@ -121,7 +116,7 @@ println("Learned noise parameter, σ_1:") learned_σ_1 = sqrt(emulator_1.machine_learning_tool.models[1].kernel.kright.σ2) println("σ_1 = $learned_σ_1") # Check if the learned noise is approximately 1 -@assert(isapprox(learned_σ_1, 1.0; atol=0.1)) +@assert(isapprox(learned_σ_1, 1.0; atol = 0.1)) println("\nKernel of the GP trained to predict y2 from x=(x1, x2):") println(emulator_1.machine_learning_tool.models[2].kernel) @@ -129,24 +124,19 @@ println("Learned noise parameter, σ_2:") learned_σ_2 = sqrt(emulator_1.machine_learning_tool.models[2].kernel.kright.σ2) println("σ_2 = $learned_σ_2") # Check if the learned noise is approximately 1 -@assert(isapprox(learned_σ_2, 1.0; atol=0.1)) +@assert(isapprox(learned_σ_2, 1.0; atol = 0.1)) println("------------------------------------------------------------------\n") # For comparison: When noise_learn is set to false, the observational noise # is set to 1.0 and is not learned/optimized during the training. But thanks # to the SVD, 1.0 is the correct value to use. gauss_proc_2 = GaussianProcess( - gppackage; - kernel=nothing, # use default squared exponential kernel - prediction_type=prediction_type, - noise_learn=false -) -emulator_2 = Emulator( - gauss_proc_2, - iopairs, - obs_noise_cov=Σ, - normalize_inputs=true + gppackage; + kernel = nothing, # use default squared exponential kernel + prediction_type = prediction_type, + noise_learn = false, ) +emulator_2 = Emulator(gauss_proc_2, iopairs, obs_noise_cov = Σ, normalize_inputs = true) println("build GP with ", n, " training points") #optimize the hyperparameters to best fit the GP. optimize_hyperparameters!(emulator_2) @@ -159,7 +149,7 @@ println("-----------") println("\nKernel of the GP trained to predict y1 from x=(x1, x2):") # Note: In contrast to the kernels of the gpobj1 models, these ones do not # have a white noise ("Noise") kernel component -println(emulator_2.machine_learning_tool.models[1].kernel) +println(emulator_2.machine_learning_tool.models[1].kernel) # logNoise is given as log(sqrt(noise)) obs_noise_1 = exp(emulator_2.machine_learning_tool.models[1].logNoise.value^2) println("Observational noise: $obs_noise_1") diff --git a/examples/Emulator/GaussianProcess/plot_GP.jl b/examples/Emulator/GaussianProcess/plot_GP.jl index cbf5da877..6c6127f8a 100644 --- a/examples/Emulator/GaussianProcess/plot_GP.jl +++ b/examples/Emulator/GaussianProcess/plot_GP.jl @@ -13,11 +13,10 @@ using CalibrateEmulateSample.DataContainers plot_flag = true if plot_flag using Plots - gr(size=(1500, 700)) + gr(size = (1500, 700)) Plots.scalefontsizes(1.3) font = Plots.font("Helvetica", 18) - fontdict = Dict(:guidefont=>font, :xtickfont=>font, :ytickfont=>font, - :legendfont=>font) + fontdict = Dict(:guidefont => font, :xtickfont => font, :ytickfont => font, :legendfont => font) end @@ -53,7 +52,7 @@ end # # ############################################################################### -function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where T +function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where {T} m, n = length(vy), length(vx) gx = reshape(repeat(vx, inner = m, outer = 1), m, n) gy = reshape(repeat(vy, inner = 1, outer = n), m, n) @@ -72,7 +71,7 @@ end #create the machine learning tools: Gaussian Process gppackage = GPJL() pred_type = YType() -gaussian_process = GaussianProcess(gppackage,noise_learn=true) +gaussian_process = GaussianProcess(gppackage, noise_learn = true) # Generate training data (x-y pairs, where x ∈ ℝ ᵖ, y ∈ ℝ ᵈ) # x = [x1, x2]: inputs/predictors/features/parameters @@ -88,48 +87,80 @@ X = 2.0 * π * rand(p, n) # G(x1, x2) g1x = sin.(X[1, :]) .+ cos.(X[2, :]) g2x = sin.(X[1, :]) .- cos.(X[2, :]) -gx = zeros(2,n) -gx[1,:] = g1x -gx[2,:] = g2x - +gx = zeros(2, n) +gx[1, :] = g1x +gx[2, :] = g2x + # Add noise η -μ = zeros(d) +μ = zeros(d) Σ = 0.1 * [[0.8, 0.0] [0.0, 0.5]] # d x d -noise_samples = rand(MvNormal(μ, Σ), n) +noise_samples = rand(MvNormal(μ, Σ), n) # y = G(x) + η Y = gx .+ noise_samples #plot training data with and without noise if plot_flag - p1 = plot(X[1,:], X[2,:], g1x, st=:surface, camera=(30, 60), c=:cividis, - xlabel="x1", ylabel="x2", - zguidefontrotation=90) - + p1 = plot( + X[1, :], + X[2, :], + g1x, + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zguidefontrotation = 90, + ) + figpath = joinpath(output_directory, "GP_test_observed_y1nonoise.png") savefig(figpath) - - p2 = plot(X[1,:], X[2,:], g2x, st=:surface, camera=(30, 60), c=:cividis, - xlabel="x1", ylabel="x2", - zguidefontrotation=90) + + p2 = plot( + X[1, :], + X[2, :], + g2x, + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zguidefontrotation = 90, + ) figpath = joinpath(output_directory, "GP_test_observed_y2nonoise.png") savefig(figpath) - - p3 = plot(X[1,:], X[2,:], Y[1,:], st=:surface, camera=(30, 60), c=:cividis, - xlabel="x1", ylabel="x2", - zguidefontrotation=90) + + p3 = plot( + X[1, :], + X[2, :], + Y[1, :], + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zguidefontrotation = 90, + ) figpath = joinpath(output_directory, "GP_test_observed_y1.png") savefig(figpath) - - p4 = plot(X[1,:], X[2,:], Y[2,:], st=:surface, camera=(30, 60), c=:cividis, - xlabel="x1", ylabel="x2", - zguidefontrotation=90) + + p4 = plot( + X[1, :], + X[2, :], + Y[2, :], + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zguidefontrotation = 90, + ) figpath = joinpath(output_directory, "GP_test_observed_y2.png") savefig(figpath) end -iopairs = PairedDataContainer(X,Y,data_are_columns=true) +iopairs = PairedDataContainer(X, Y, data_are_columns = true) @assert get_inputs(iopairs) == X @assert get_outputs(iopairs) == Y @@ -140,11 +171,7 @@ iopairs = PairedDataContainer(X,Y,data_are_columns=true) # exponential kernel. # Setting noise_learn=true leads to the addition of white noise to the # kernel -emulator = Emulator( - gaussian_process, - iopairs, - obs_noise_cov=Σ, - normalize_inputs=false) +emulator = Emulator(gaussian_process, iopairs, obs_noise_cov = Σ, normalize_inputs = false) println("build GP with ", n, " training points") #optimize the hyperparameters to best fit the GP. @@ -157,11 +184,11 @@ println("GP trained") # Plot mean and variance of the predicted observables y1 and y2 # For this, we generate test points on a x1-x2 grid. n_pts = 50 -x1 = range(0.0, stop=2*π, length=n_pts) -x2 = range(0.0, stop=2*π, length=n_pts) +x1 = range(0.0, stop = 2 * π, length = n_pts) +x2 = range(0.0, stop = 2 * π, length = n_pts) X1, X2 = meshgrid(x1, x2) # Input for predict has to be of size N_samples x input_dim -inputs = permutedims(hcat(X1[:], X2[:]),(2,1)) +inputs = permutedims(hcat(X1[:], X2[:]), (2, 1)) #gp_mean, gp_cov = GaussianProcessEmulator.predict(gpobj, # inputs, @@ -170,52 +197,87 @@ inputs = permutedims(hcat(X1[:], X2[:]),(2,1)) # covariance matrices, not just the variance -- gp_cov is a vector # of covariance matrices) -gp_mean, gp_cov = predict(emulator, inputs, transform_to_real=true) -println("end predictions at ", n_pts*n_pts, " points") +gp_mean, gp_cov = predict(emulator, inputs, transform_to_real = true) +println("end predictions at ", n_pts * n_pts, " points") #plot predictions for y_i in 1:d - gp_var_temp = [diag(gp_cov[j]) for j in 1:length(gp_cov)] # (40000,) - gp_var = permutedims(vcat([x' for x in gp_var_temp]...),(2,1)) # 2 x 40000 - + gp_var = permutedims(vcat([x' for x in gp_var_temp]...), (2, 1)) # 2 x 40000 + mean_grid = reshape(gp_mean[y_i, :], n_pts, n_pts) # 2 x 40000 if plot_flag - p5 = plot(x1, x2, mean_grid, st=:surface, camera=(30, 60), c=:cividis, - xlabel="x1", ylabel="x2", zlabel="mean of y"*string(y_i), - zguidefontrotation=90) + p5 = plot( + x1, + x2, + mean_grid, + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zlabel = "mean of y" * string(y_i), + zguidefontrotation = 90, + ) end var_grid = reshape(gp_var[y_i, :], n_pts, n_pts) if plot_flag - p6 = plot(x1, x2, var_grid, st=:surface, camera=(30, 60), c=:cividis, - xlabel="x1", ylabel="x2", zlabel="var of y"*string(y_i), - zguidefontrotation=90) - - plot(p5, p6, layout=(1, 2), legend=false) - - savefig(joinpath(output_directory, "GP_test_y"*string(y_i)*"_predictions.png")) + p6 = plot( + x1, + x2, + var_grid, + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zlabel = "var of y" * string(y_i), + zguidefontrotation = 90, + ) + + plot(p5, p6, layout = (1, 2), legend = false) + + savefig(joinpath(output_directory, "GP_test_y" * string(y_i) * "_predictions.png")) end end # Plot the true components of G(x1, x2) -g1_true = sin.(inputs[1,:]) .+ cos.(inputs[2, :]) +g1_true = sin.(inputs[1, :]) .+ cos.(inputs[2, :]) g1_true_grid = reshape(g1_true, n_pts, n_pts) if plot_flag - p7 = plot(x1, x2, g1_true_grid, st=:surface, camera=(30, 60), c=:cividis, - xlabel="x1", ylabel="x2", zlabel="sin(x1) + cos(x2)", - zguidefontrotation=90) + p7 = plot( + x1, + x2, + g1_true_grid, + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zlabel = "sin(x1) + cos(x2)", + zguidefontrotation = 90, + ) savefig(joinpath(output_directory, "GP_test_true_g1.png")) end g2_true = sin.(inputs[1, :]) .- cos.(inputs[2, :]) g2_true_grid = reshape(g2_true, n_pts, n_pts) if plot_flag - p8 = plot(x1, x2, g2_true_grid, st=:surface, camera=(30, 60), c=:cividis, - xlabel="x1", ylabel="x2", zlabel="sin(x1) - cos(x2)", - zguidefontrotation=90) + p8 = plot( + x1, + x2, + g2_true_grid, + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zlabel = "sin(x1) - cos(x2)", + zguidefontrotation = 90, + ) g_true_grids = [g1_true_grid, g2_true_grid] - + savefig(joinpath(output_directory, "GP_test_true_g2.png")) end @@ -225,21 +287,29 @@ for y_i in 1:d # Reshape gp_cov to size N_samples x output_dim gp_var_temp = [diag(gp_cov[j]) for j in 1:length(gp_cov)] # (40000,) - gp_var = permutedims(vcat([x' for x in gp_var_temp]...),(2,1)) # 40000 x 2 + gp_var = permutedims(vcat([x' for x in gp_var_temp]...), (2, 1)) # 40000 x 2 - mean_grid = reshape(gp_mean[y_i,:], n_pts, n_pts) + mean_grid = reshape(gp_mean[y_i, :], n_pts, n_pts) var_grid = reshape(gp_var[y_i, :], n_pts, n_pts) # Compute and plot 1/variance * (truth - prediction)^2 if plot_flag - zlabel = "1/var * (true_y"*string(y_i)*" - predicted_y"*string(y_i)*")^2" - - p9 = plot(x1, x2, sqrt.(1.0 ./ var_grid .* (g_true_grids[y_i] .- mean_grid).^2), - st=:surface, camera=(30, 60), c=:magma, zlabel=zlabel, - xlabel="x1", ylabel="x2", - zguidefontrotation=90) - - savefig(joinpath(output_directory, "GP_test_y"*string(y_i)*"_difference_truth_prediction.png")) + zlabel = "1/var * (true_y" * string(y_i) * " - predicted_y" * string(y_i) * ")^2" + + p9 = plot( + x1, + x2, + sqrt.(1.0 ./ var_grid .* (g_true_grids[y_i] .- mean_grid) .^ 2), + st = :surface, + camera = (30, 60), + c = :magma, + zlabel = zlabel, + xlabel = "x1", + ylabel = "x2", + zguidefontrotation = 90, + ) + + savefig(joinpath(output_directory, "GP_test_y" * string(y_i) * "_difference_truth_prediction.png")) end end diff --git a/examples/Lorenz/GModel.jl b/examples/Lorenz/GModel.jl index c7b5d4491..f76edcbd1 100644 --- a/examples/Lorenz/GModel.jl +++ b/examples/Lorenz/GModel.jl @@ -14,118 +14,110 @@ export run_G_ensemble export lorenz_forward """ - GSettings{FT<:AbstractFloat, KT, D} +$(DocStringExtensions.TYPEDEF) -Structure to hold all information to run the forward model G +Structure to hold all information to run the forward model *G*. # Fields -$(DocStringExtensions.FIELDS) +$(DocStringExtensions.TYPEDFIELDS) """ -struct GSettings{FT<:AbstractFloat, KT, D} - "model output" +struct GSettings{FT <: AbstractFloat, KT, D} + "Model output." out::Array{FT, 1} - "time period over which to run the model, e.g., `(0, 1)`" + "Time period over which to run the model, e.g., `(0, 1)`." tspan::Tuple{FT, FT} - "x domain" + "X domain." x::Array{FT, 1} - "y domain" + "Y domain." y::Array{FT, 1} end -struct LSettings - # Stationary or transient dynamics - dynamics::Int32 - # G model statistics type - stats_type::Int32 - # Integration time start - t_start::Float64 - # Data collection length - T::Float64 - # Integration length - Ts::Float64 - # Duration of polynomial fit, number of samples of Ts to fit over - Tfit::Float64 - # Initial perturbation - Fp::Array{Float64} - # Number of longitude steps - N::Int32 - # Timestep - dt::Float64 - # Simulation end time - tend::Float64 - # For stats_type=2, number of frequencies to consider - kmax::Int32 - # Should ω be learned using CES? - ω_fixed::Bool - # Truth ω - ω_true::Float64 +struct LSettings + # Stationary or transient dynamics + dynamics::Int32 + # G model statistics type + stats_type::Int32 + # Integration time start + t_start::Float64 + # Data collection length + T::Float64 + # Integration length + Ts::Float64 + # Duration of polynomial fit, number of samples of Ts to fit over + Tfit::Float64 + # Initial perturbation + Fp::Array{Float64} + # Number of longitude steps + N::Int32 + # Timestep + dt::Float64 + # Simulation end time + tend::Float64 + # For stats_type=2, number of frequencies to consider + kmax::Int32 + # Should ω be learned using CES? + ω_fixed::Bool + # Truth ω + ω_true::Float64 end struct LParams - # Mean forcing - F::Float64 - # Forcing frequency for transient terms - ω::Float64 - # Forcing amplitude for transient terms - A::Float64 + # Mean forcing + F::Float64 + # Forcing frequency for transient terms + ω::Float64 + # Forcing amplitude for transient terms + A::Float64 end """ - run_G_ensemble(params::Array{FT, 2}, - settings::GSettings{FT}, - update_params, - moment, - get_src; - rng_seed=42) where {FT<:AbstractFloat} +$(DocStringExtensions.TYPEDSIGNATURES) -Run the forward model G for an array of parameters by iteratively -calling run_G for each of the N_ensemble parameter values. -Return g_ens, an array of size N_data x N_ensemble, where -g_ens[:,j] = G(params[:,j]) +Run the forward model *G* for an array of parameters by iteratively +calling `run_G` for each of the *N\\_ensemble* parameter values. - - `params` - array of size N_parameters x N_ensemble containing the - parameters for which G will be run - - `settings` - a GSetttings struct +- `params` - array of size (*N\\_ensemble* × *N\\_parameters*) containing the parameters for + which *G* will be run. +- `settings` - a [GSetttings](@ref) struct. +Returns `g_ens`, an array of size (*N\\_ensemble* × *N\\_data*), where g_ens[j,:] = G(params[j,:]). """ -function run_G_ensemble(params::Array{FT, 2}, - settings::LSettings, - rng_seed=42) where {FT<:AbstractFloat} +function run_G_ensemble(params::Array{FT, 2}, settings::LSettings, rng_seed = 42) where {FT <: AbstractFloat} # Initilize ensemble N_ens = size(params, 2) # params is N_params x N_ens if settings.stats_type == 1 - nd = 2 + nd = 2 elseif settings.stats_type == 2 - nd = 1 + (2*settings.kmax) + nd = 1 + (2 * settings.kmax) elseif settings.stats_type == 3 - nd = 3 + nd = 3 elseif settings.stats_type == 4 - nd = Int64(2*(settings.T/settings.Ts/settings.Tfit)) + nd = Int64(2 * (settings.T / settings.Ts / settings.Tfit)) elseif settings.stats_type == 5 - nd = Int64(settings.T/settings.Ts/settings.Tfit) + nd = Int64(settings.T / settings.Ts / settings.Tfit) end g_ens = zeros(nd, N_ens) # Lorenz parameters #Fp = rand(Normal(0.0, 0.01), N); # Initial perturbation F = params[1, :] # Forcing - if settings.dynamics==2 + if settings.dynamics == 2 A = params[2, :] # Transience amplitude else - A = F .* 0. # Set to zero, it is not used + A = F .* 0.0 # Set to zero, it is not used end - if settings.ω_fixed==false; - ω = params[3, :] # Transience frequency + if settings.ω_fixed == false + ω = params[3, :] # Transience frequency end Random.seed!(rng_seed) for i in 1:N_ens # run the model with the current parameters, i.e., map θ to G(θ) - if settings.ω_fixed==false - lorenz_params = GModel.LParams(F[i], ω[i], A[i]) - elseif settings.ω_fixed==true - lorenz_params = GModel.LParams(F[i], settings.ω_true, A[i]) - end - g_ens[:, i] = lorenz_forward(settings, lorenz_params) + if settings.ω_fixed == false + lorenz_params = GModel.LParams(F[i], ω[i], A[i]) + elseif settings.ω_fixed == true + lorenz_params = GModel.LParams(F[i], settings.ω_true, A[i]) + end + g_ens[:, i] = lorenz_forward(settings, lorenz_params) end return g_ens @@ -136,12 +128,12 @@ end # F: scalar forcing Float64(1), Fp: initial perturbation Float64(N) # N: number of longitude steps Int64(1), dt: time step Float64(1) # tend: End of simulation Float64(1), nstep: -function lorenz_forward(settings::LSettings, params::LParams) - # run the Lorenz simulation - xn, t = lorenz_solve(settings, params) - # Get statistics - gt = stats(settings, xn, t) - return gt +function lorenz_forward(settings::LSettings, params::LParams) + # run the Lorenz simulation + xn, t = lorenz_solve(settings, params) + # Get statistics + gt = stats(settings, xn, t) + return gt end @@ -150,202 +142,223 @@ end # N: number of longitude steps Int64(1), dt: time step Float64(1) # tend: End of simulation Float64(1), nstep: function lorenz_solve(settings::LSettings, params::LParams) - # Initialize - nstep = Int32(ceil(settings.tend/settings.dt)); - xn = zeros(Float64, settings.N, nstep) - t = zeros(Float64, nstep) - # Initial perturbation - X = fill(Float64(0.), settings.N) - X = X + settings.Fp - # March forward in time - for j in 1:nstep - t[j] = settings.dt*j - if settings.dynamics==1 - X = RK4(X,settings.dt,settings.N,params.F) - elseif settings.dynamics==2 - F_local = params.F + params.A*sin(params.ω*t[j]) - X = RK4(X,settings.dt,settings.N,F_local) - end - xn[:,j] = X - end - # Output - return xn, t + # Initialize + nstep = Int32(ceil(settings.tend / settings.dt)) + xn = zeros(Float64, settings.N, nstep) + t = zeros(Float64, nstep) + # Initial perturbation + X = fill(Float64(0.0), settings.N) + X = X + settings.Fp + # March forward in time + for j in 1:nstep + t[j] = settings.dt * j + if settings.dynamics == 1 + X = RK4(X, settings.dt, settings.N, params.F) + elseif settings.dynamics == 2 + F_local = params.F + params.A * sin(params.ω * t[j]) + X = RK4(X, settings.dt, settings.N, F_local) + end + xn[:, j] = X + end + # Output + return xn, t -end +end # Lorenz 96 system # f = dx/dt # Inputs: x: state, N: longitude steps, F: forcing -function f(x,N,F) - f = zeros(Float64, N) - # Loop over N positions - for i in 3:N-1 - f[i] = -x[i-2]*x[i-1] + x[i-1]*x[i+1] - x[i] + F - end - # Periodic boundary conditions - f[1] = -x[N-1]*x[N] + x[N]*x[2] - x[1] + F - f[2] = -x[N]*x[1] + x[1]*x[3] - x[2] + F - f[N] = -x[N-2]*x[N-1] + x[N-1]*x[1] - x[N] + F - # Output - return f +function f(x, N, F) + f = zeros(Float64, N) + # Loop over N positions + for i in 3:(N - 1) + f[i] = -x[i - 2] * x[i - 1] + x[i - 1] * x[i + 1] - x[i] + F + end + # Periodic boundary conditions + f[1] = -x[N - 1] * x[N] + x[N] * x[2] - x[1] + F + f[2] = -x[N] * x[1] + x[1] * x[3] - x[2] + F + f[N] = -x[N - 2] * x[N - 1] + x[N - 1] * x[1] - x[N] + F + # Output + return f end # RK4 solve function RK4(xold, dt, N, F) - # Predictor steps - k1 = f(xold, N, F) - k2 = f(xold + k1*dt/2., N, F) - k3 = f(xold + k2*dt/2., N, F) - k4 = f(xold + k3*dt, N, F) - # Step - xnew = xold + (dt/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4) - # Output - return xnew + # Predictor steps + k1 = f(xold, N, F) + k2 = f(xold + k1 * dt / 2.0, N, F) + k3 = f(xold + k2 * dt / 2.0, N, F) + k4 = f(xold + k3 * dt, N, F) + # Step + xnew = xold + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4) + # Output + return xnew end -function regression(X,y) - beta_tilde = [ ones(length(X),1) X ] \ y; - v = beta_tilde[1]; beta = beta_tilde[2]; - return beta, v +function regression(X, y) + beta_tilde = [ones(length(X), 1) X] \ y + v = beta_tilde[1] + beta = beta_tilde[2] + return beta, v end -function spectra(signal, Ts, t) - # Init - Nt = length(t) - if mod(Nt,2)!=0 - t=t[1:Nt-1]; signal = signal[1:Nt-1]; Nt = Nt-1; - end - # Remove mean - u = signal .- mean(signal,dims=1) - # FFT - F = fft(u) |> fftshift - freqs = fftfreq(length(t), 1.0/Ts) |> fftshift - #k = gen_k(Nt,Ts) |> fftshift - T = (Nt-1)*Ts; - A = Ts^2 / (2*pi*T); - uhat = A*abs.(F.^2); - # Output - mid = Int32(Nt/2)+1; - f = freqs[mid:Nt] - Fp = uhat[mid:Nt] - return f, Fp +function spectra(signal, Ts, t) + # Init + Nt = length(t) + if mod(Nt, 2) != 0 + t = t[1:(Nt - 1)] + signal = signal[1:(Nt - 1)] + Nt = Nt - 1 + end + # Remove mean + u = signal .- mean(signal, dims = 1) + # FFT + F = fft(u) |> fftshift + freqs = fftfreq(length(t), 1.0 / Ts) |> fftshift + #k = gen_k(Nt,Ts) |> fftshift + T = (Nt - 1) * Ts + A = Ts^2 / (2 * pi * T) + uhat = A * abs.(F .^ 2) + # Output + mid = Int32(Nt / 2) + 1 + f = freqs[mid:Nt] + Fp = uhat[mid:Nt] + return f, Fp end ############################### ## Extract statistics ############################### -function stats(settings, xn, t) -# Define averaging indices range - indices = findall(x -> (x>settings.t_start) - && (x (x > settings.t_start) && (x < settings.t_start + settings.T), t) + # Define statistics of interest + if settings.stats_type == 1 # Mean + # Average in time and over longitude + gtm = mean(mean(xn[:, indices], dims = 2), dims = 1) + # Variance + gtc = mean(mean((xn[:, indices] .- gtm) .^ 2, dims = 2), dims = 1) + # Combine statistics + gt = vcat(gtm, gtc) + elseif settings.stats_type == 2 #What about an integral under parts of the spectrum + # Average in time and over longitude + gtm = mean(mean(xn[:, indices], dims = 2), dims = 1) + # Power spectra + # Get size + f, Fp = spectra(xn[1, indices]', settings.dt, t[indices]) + # Compute temporal spectrum for each longitude + Fp = zeros(size(Fp, 1), settings.N) + for i in 1:(settings.N) + f, Fp[:, i] = spectra(xn[i, indices]', settings.dt, t[indices]) + end + # Average spectra over periodic directions + Fp = dropdims(mean(Fp, dims = 2), dims = 2) + ys = partialsortperm(Fp, 1:(settings.kmax); rev = true) + gt = vcat(gtm, Fp[ys]..., f[ys]...) + elseif settings.stats_type == 3 # Structure function + # Average in time and over longitude + gtm = mean(mean(xn[:, indices], dims = 2), dims = 1) + # Maximum + mxval, mxind = findmax(xn[:, indices]; dims = 2) + mxval_out = mean(mxval, dims = 1) + # Power spectra + # Get size + f, Fp = spectra(xn[1, indices]', settings.dt, t[indices]) + # Compute temporal spectrum for each longitude + Fp = zeros(size(Fp, 1), settings.N) + for i in 1:(settings.N) + f, Fp[:, i] = spectra(xn[i, indices]', settings.dt, t[indices]) + end + # Average spectra over periodic directions + Fp = dropdims(mean(Fp, dims = 2), dims = 2) + ys = partialsortperm(Fp, 1; rev = true) + # Period + T = 1.0 / f[ys] + mxval, r = findmin(abs.(t .- T); dims = 1) + r = r[1] + # Structure function + xp = xn .- gtm + st = zeros(settings.N) + for i in 1:(settings.N) + st[i] = mean( + ( + xp[1, indices[1]:(indices[length(indices)] - r)] .- + xp[1, (indices[1] + r):indices[length(indices)]] + ) .^ 2, + ) + end + # Combine + gt = vcat(gtm, mean(st)..., mxval_out...) + elseif settings.stats_type == 4 # Variance sub-samples, linear fit over window + T = settings.T + # Calculate variance on subsamples + # Variances over sliding windows + Ts = settings.Ts + N = Int64(T / Ts) # number of sliding windows + nt = Int64(Ts / settings.dt) - 1 # length of sliding windows + var_slide = zeros(N) + t_slide = zeros(nt, N) + for i in 1:N + var_slide[i] = var(vcat(xn[:, indices[((i - 1) * nt + 1):(i * nt)]]...)) + t_slide[:, i] = t[((i - 1) * nt + 1):(i * nt)] + end + # Polynomial fit over a batch of sliding windows + Tfit = Int64(settings.Tfit) + NT = Int64(N / Tfit) + a1 = zeros(NT) + a2 = zeros(NT) + t_start = zeros(NT) + y = zeros(Tfit, NT) + ty = zeros(Tfit, NT) + ym = zeros(Tfit, NT) + for i in 1:NT + ind_low = (i - 1) * Tfit + 1 + ind_high = i * Tfit + t_start[i] = t_slide[1, ind_low] + ty[:, i] = t_slide[1, ind_low:ind_high] .- t_start[i] + a1[i], a2[i] = regression(ty[:, i], var_slide[ind_low:ind_high]) + y[:, i] = a1[i] .* ty[:, i] .+ a2[i] + ym[i] = mean(var_slide[ind_low:ind_high]) + end + # Combine + gt = vcat(a1..., a2...) + elseif settings.stats_type == 5 # Variance sub-samples, mean over window + T = settings.T + # Calculate variance on subsamples + # Variances over sliding windows + Ts = settings.Ts + N = Int64(T / Ts) # number of sliding windows + nt = Int64(Ts / settings.dt) - 1 # length of sliding windows + var_slide = zeros(N) + t_slide = zeros(nt, N) + for i in 1:N + var_slide[i] = var(vcat(xn[:, indices[((i - 1) * nt + 1):(i * nt)]]...)) + t_slide[:, i] = t[((i - 1) * nt + 1):(i * nt)] + end + # Polynomial fit over a batch of sliding windows + Tfit = Int64(settings.Tfit) + NT = Int64(N / Tfit) + a1 = zeros(NT) + a2 = zeros(NT) + t_start = zeros(NT) + y = zeros(Tfit, NT) + ty = zeros(Tfit, NT) + ym = zeros(NT) + for i in 1:NT + ind_low = (i - 1) * Tfit + 1 + ind_high = i * Tfit + t_start[i] = t_slide[1, ind_low] + ty[:, i] = t_slide[1, ind_low:ind_high] .- t_start[i] + a1[i], a2[i] = regression(ty[:, i], var_slide[ind_low:ind_high]) + y[:, i] = a1[i] .* ty[:, i] .+ a2[i] + ym[i] = mean(var_slide[ind_low:ind_high]) + end + # Combine + gt = vcat(ym...) + else + ArgumentError("Setting " * string(settings.stats_type) * " not implemented.") + end + return gt end diff --git a/examples/Lorenz/Lorenz_example.jl b/examples/Lorenz/Lorenz_example.jl index f31f7699f..1243e9bae 100644 --- a/examples/Lorenz/Lorenz_example.jl +++ b/examples/Lorenz/Lorenz_example.jl @@ -29,14 +29,14 @@ rng_seed = 2413798 Random.seed!(rng_seed) -function get_standardizing_factors(data::Array{FT,2}) where {FT} +function get_standardizing_factors(data::Array{FT, 2}) where {FT} # Input: data size: N_data x N_ensembles # Ensemble median of the data - norm_factor = median(data,dims=2) + norm_factor = median(data, dims = 2) return norm_factor end -function get_standardizing_factors(data::Array{FT,1}) where {FT} +function get_standardizing_factors(data::Array{FT, 1}) where {FT} # Input: data size: N_data*N_ensembles (splatted) # Ensemble median of the data norm_factor = median(data) @@ -58,27 +58,27 @@ end # Governing settings # Characteristic time scale -τc = 5. # days, prescribed by the L96 problem +τc = 5.0 # days, prescribed by the L96 problem # Stationary or transient dynamics dynamics = 2 # Transient is 2 # Statistics integration length # This has to be less than 360 and 360 must be divisible by Ts_days -Ts_days = 30. # Integration length in days +Ts_days = 30.0 # Integration length in days # Stats type, which statistics to construct from the L96 system # 4 is a linear fit over a batch of length Ts_days # 5 is the mean over a batch of length Ts_days -stats_type = 5 +stats_type = 5 ### ### Define the (true) parameters ### # Define the parameters that we want to learn -F_true = 8. # Mean F +F_true = 8.0 # Mean F A_true = 2.5 # Transient F amplitude -ω_true = 2. * π / (360. / τc) # Frequency of the transient F +ω_true = 2.0 * π / (360.0 / τc) # Frequency of the transient F -if dynamics==2 +if dynamics == 2 params_true = [F_true, A_true] param_names = ["F", "A"] else @@ -86,7 +86,7 @@ else param_names = ["F"] end n_param = length(param_names) -params_true = reshape(params_true, (n_param,1)) +params_true = reshape(params_true, (n_param, 1)) println(n_param) println(params_true) @@ -99,30 +99,30 @@ println(params_true) log_normal = false # THIS ISN't CURRENTLY IMPLEMENTED function logmean_and_logstd(μ, σ) - σ_log = sqrt(log(1.0 + σ^2/μ^2)) - μ_log = log(μ / (sqrt(1.0 + σ^2/μ^2))) + σ_log = sqrt(log(1.0 + σ^2 / μ^2)) + μ_log = log(μ / (sqrt(1.0 + σ^2 / μ^2))) return μ_log, σ_log end #logmean_F, logstd_F = logmean_and_logstd(F_true, 5) #logmean_A, logstd_A = logmean_and_logstd(A_true, 0.2*A_true) if dynamics == 2 - #prior_means = [F_true+0.5, A_true+0.5] - prior_means = [F_true, A_true] - prior_stds = [3.0, 0.5*A_true] - d1 = Parameterized(Normal(prior_means[1], prior_stds[1])) - d2 = Parameterized(Normal(prior_means[2], prior_stds[2])) - prior_distns = [d1, d2] - c1 = no_constraint() - c2 = no_constraint() - constraints = [[c1], [c2]] - prior_names = param_names - priors = ParameterDistribution(prior_distns, constraints, prior_names) + #prior_means = [F_true+0.5, A_true+0.5] + prior_means = [F_true, A_true] + prior_stds = [3.0, 0.5 * A_true] + d1 = Parameterized(Normal(prior_means[1], prior_stds[1])) + d2 = Parameterized(Normal(prior_means[2], prior_stds[2])) + prior_distns = [d1, d2] + c1 = no_constraint() + c2 = no_constraint() + constraints = [[c1], [c2]] + prior_names = param_names + priors = ParameterDistribution(prior_distns, constraints, prior_names) else - prior_distns = [Parameterized(Normal(F_true, 1))] - constraints = [[no_constraint()]] - prior_names = ["F"] - priors = ParameterDistribution(prior_distns, constraints, prior_names) + prior_distns = [Parameterized(Normal(F_true, 1))] + constraints = [[no_constraint()]] + prior_names = ["F"] + priors = ParameterDistribution(prior_distns, constraints, prior_names) end @@ -139,18 +139,18 @@ data_names = ["y0", "y1"] # Lorenz 96 model parameters # Behavior changes depending on the size of F, N N = 36 -dt = 1/64. +dt = 1 / 64.0 # Start of integration -t_start = 800. +t_start = 800.0 # Data collection length #if dynamics==1 # T = 2. #else # T = 360. / τc #end -T = 360. / τc +T = 360.0 / τc # Batch length -Ts = 5. / τc # Nondimensionalize by L96 timescale +Ts = 5.0 / τc # Nondimensionalize by L96 timescale # Integration length Tfit = Ts_days / τc # Initial perturbation @@ -163,9 +163,8 @@ var_prescribe = false # Settings # Constructs an LSettings structure, see GModel.jl for the descriptions -lorenz_settings = GModel.LSettings(dynamics, stats_type, t_start, - T, Ts, Tfit, Fp, N, dt, t_start+T, kmax, - ω_fixed, ω_true); +lorenz_settings = + GModel.LSettings(dynamics, stats_type, t_start, T, Ts, Tfit, Fp, N, dt, t_start + T, kmax, ω_fixed, ω_true); lorenz_params = GModel.LParams(F_true, ω_true, A_true) ### @@ -179,34 +178,46 @@ lorenz_params = GModel.LParams(F_true, ω_true, A_true) # Output: gt: [N_data, N_ens] # Dropdims of the output since the forward model is only being run with N_ens=1 # corresponding to the truth construction -gt = dropdims(GModel.run_G_ensemble(params_true, lorenz_settings), dims=2) +gt = dropdims(GModel.run_G_ensemble(params_true, lorenz_settings), dims = 2) # Compute internal variability covariance -n_samples = 50 -if var_prescribe==true +n_samples = 50 +if var_prescribe == true n_samples = 100 - yt = zeros(length(gt),n_samples) + yt = zeros(length(gt), n_samples) noise_level = 0.05 Γy = noise_level * convert(Array, Diagonal(gt)) μ = zeros(length(gt)) # Add noise for i in 1:n_samples - yt[:,i] = gt .+ rand(MvNormal(μ, Γy)) + yt[:, i] = gt .+ rand(MvNormal(μ, Γy)) end else println("Using truth values to compute covariance") yt = zeros(length(gt), n_samples) for i in 1:n_samples - lorenz_settings_local = GModel.LSettings(dynamics, stats_type, t_start+T*(i-1), - T, Ts, Tfit, Fp, N, dt, - t_start+T*(i-1)+T, kmax, ω_fixed, ω_true); - yt[:, i] = GModel.run_G_ensemble(params_true, lorenz_settings_local) + lorenz_settings_local = GModel.LSettings( + dynamics, + stats_type, + t_start + T * (i - 1), + T, + Ts, + Tfit, + Fp, + N, + dt, + t_start + T * (i - 1) + T, + kmax, + ω_fixed, + ω_true, + ) + yt[:, i] = GModel.run_G_ensemble(params_true, lorenz_settings_local) end # Variance of truth data #Γy = convert(Array, Diagonal(dropdims(mean((yt.-mean(yt,dims=1)).^2,dims=1),dims=1))) # Covariance of truth data - Γy = cov(yt, dims=2) - + Γy = cov(yt, dims = 2) + println(Γy) println(size(Γy), " ", rank(Γy)) end @@ -235,29 +246,26 @@ exp_transform(a::AbstractArray) = exp.(a) N_ens = 50 # number of ensemble members N_iter = 5 # number of EKI iterations # initial parameters: N_params x N_ens -initial_params = construct_initial_ensemble(priors, N_ens; rng_seed=rng_seed) +initial_params = construct_initial_ensemble(priors, N_ens; rng_seed = rng_seed) -ekiobj = EnsembleKalmanProcesses.EnsembleKalmanProcess(initial_params, - truth_sample, - truth.obs_noise_cov, - Inversion()) +ekiobj = EnsembleKalmanProcesses.EnsembleKalmanProcess(initial_params, truth_sample, truth.obs_noise_cov, Inversion()) # EKI iterations println("EKP inversion error:") -err = zeros(N_iter) -err_params = zeros(N_iter) +err = zeros(N_iter) +err_params = zeros(N_iter) for i in 1:N_iter - if log_normal==false + if log_normal == false params_i = get_u_final(ekiobj) else params_i = exp_transform(get_u_final(ekiobj)) end g_ens = GModel.run_G_ensemble(params_i, lorenz_settings_G) - EnsembleKalmanProcesses.update_ensemble!(ekiobj, g_ens) - err[i] = get_error(ekiobj)[end] - err_params[i] = mean((params_true - mean(params_i,dims=2)).^2) - println("Iteration: "*string(i)*", Error (data): "*string(err[i])) - println("Iteration: "*string(i)*", Error (params): "*string(err_params[i])) + EnsembleKalmanProcesses.update_ensemble!(ekiobj, g_ens) + err[i] = get_error(ekiobj)[end] + err_params[i] = mean((params_true - mean(params_i, dims = 2)) .^ 2) + println("Iteration: " * string(i) * ", Error (data): " * string(err[i])) + println("Iteration: " * string(i) * ", Error (params): " * string(err_params[i])) end # EKI results: Has the ensemble collapsed toward the truth? @@ -265,10 +273,10 @@ println("True parameters: ") println(params_true) println("\nEKI results:") -if log_normal==false - println(mean(get_u_final(ekiobj), dims=2)) +if log_normal == false + println(mean(get_u_final(ekiobj), dims = 2)) else - println(mean(exp_transform(get_u_final(ekiobj)), dims=2)) + println(mean(exp_transform(get_u_final(ekiobj)), dims = 2)) end ### @@ -282,16 +290,16 @@ truncate_svd = 0.95 gppackage = Emulators.GPJL() pred_type = Emulators.YType() gauss_proc = GaussianProcess( - gppackage; - kernel=nothing, # use default squared exponential kernel - prediction_type=pred_type, - noise_learn=false + gppackage; + kernel = nothing, # use default squared exponential kernel + prediction_type = pred_type, + noise_learn = false, ) # Standardize the output data # Use median over all data since all data are the same type yt_norm = vcat(yt...) -norm_factor = get_standardizing_factors(yt_norm) +norm_factor = get_standardizing_factors(yt_norm) println(size(norm_factor)) #norm_factor = vcat(norm_factor...) norm_factor = fill(norm_factor, size(truth_sample)) @@ -305,38 +313,38 @@ normalized = true emulator = Emulator( gauss_proc, input_output_pairs; - obs_noise_cov=Γy, - normalize_inputs=normalized, - standardize_outputs=standardize, - standardize_outputs_factors=norm_factor, - truncate_svd=truncate_svd + obs_noise_cov = Γy, + normalize_inputs = normalized, + standardize_outputs = standardize, + standardize_outputs_factors = norm_factor, + truncate_svd = truncate_svd, ) optimize_hyperparameters!(emulator) # Check how well the Gaussian Process regression predicts on the # true parameters #if truncate_svd==1.0 - if log_normal==false - y_mean, y_var = Emulators.predict(emulator, params_true, transform_to_real=true) - else - y_mean, y_var = Emulators.predict(emulator, log.(params_true), transform_to_real=true) - end - - println("GP prediction on true parameters: ") - println(vec(y_mean)) - println(" GP variance") - println(diag(y_var[1],0)) - println("true data: ") - println(truth.mean) # same, regardless of norm_factor - println("GP MSE: ") - println(mean((truth.mean - vec(y_mean)).^2)) +if log_normal == false + y_mean, y_var = Emulators.predict(emulator, reshape(params_true, :, 1), transform_to_real = true) +else + y_mean, y_var = Emulators.predict(emulator, reshape(log.(params_true), :, 1), transform_to_real = true) +end + +println("GP prediction on true parameters: ") +println(vec(y_mean)) +println(" GP variance") +println(diag(y_var[1], 0)) +println("true data: ") +println(truth.mean) # same, regardless of norm_factor +println("GP MSE: ") +println(mean((truth.mean - vec(y_mean)) .^ 2)) #end ### ### Sample: Markov Chain Monte Carlo ### # initial values -u0 = vec(mean(get_inputs(input_output_pairs), dims=2)) +u0 = vec(mean(get_inputs(input_output_pairs), dims = 2)) println("initial parameters: ", u0) # MCMC parameters @@ -349,11 +357,20 @@ step = 0.1 # first guess max_iter = 2000 yt_sample = truth_sample mcmc_test = MarkovChainMonteCarlo.MCMC( - yt_sample, Γy, priors, step, u0, max_iter, mcmc_alg, burnin; - svdflag=svd_flag, standardize=standardize, truncate_svd=truncate_svd, - norm_factor=norm_factor + yt_sample, + Γy, + priors, + step, + u0, + max_iter, + mcmc_alg, + burnin; + svdflag = svd_flag, + standardize = standardize, + truncate_svd = truncate_svd, + norm_factor = norm_factor, ) -new_step = MarkovChainMonteCarlo.find_mcmc_step!(mcmc_test, emulator, max_iter=max_iter) +new_step = MarkovChainMonteCarlo.find_mcmc_step!(mcmc_test, emulator, max_iter = max_iter) # Now begin the actual MCMC println("Begin MCMC - with step size ", new_step) @@ -363,16 +380,25 @@ burnin = 2000 max_iter = 100000 mcmc = MarkovChainMonteCarlo.MCMC( - yt_sample, Γy, priors, new_step, u0, max_iter, mcmc_alg, burnin; - svdflag=svd_flag, standardize=standardize, truncate_svd=truncate_svd, - norm_factor=norm_factor + yt_sample, + Γy, + priors, + new_step, + u0, + max_iter, + mcmc_alg, + burnin; + svdflag = svd_flag, + standardize = standardize, + truncate_svd = truncate_svd, + norm_factor = norm_factor, ) MarkovChainMonteCarlo.sample_posterior!(mcmc, emulator, max_iter) println("Posterior size") println(size(mcmc.posterior)) -posterior = MarkovChainMonteCarlo.get_posterior(mcmc) +posterior = MarkovChainMonteCarlo.get_posterior(mcmc) post_mean = mean(posterior) post_cov = cov(posterior) @@ -385,54 +411,59 @@ println(det(inv(post_cov))) println(" ") # Plot the posteriors together with the priors and the true parameter values -if log_normal==false +if log_normal == false true_values = params_true else true_values = log.(params_true) end n_params = length(true_values) -y_folder = "qoi_"*string(lorenz_settings.stats_type)*"_" -save_directory = figure_save_directory*y_folder +y_folder = "qoi_" * string(lorenz_settings.stats_type) * "_" +save_directory = figure_save_directory * y_folder for idx in 1:n_params if idx == 1 param = "F" - xbounds = [F_true-1.0, F_true+1.0] - xs = collect(xbounds[1]:(xbounds[2]-xbounds[1])/100:xbounds[2]) + xbounds = [F_true - 1.0, F_true + 1.0] + xs = collect(xbounds[1]:((xbounds[2] - xbounds[1]) / 100):xbounds[2]) elseif idx == 2 param = "A" - xbounds = [A_true-1.0, A_true+1.0] - xs = collect(xbounds[1]:(xbounds[2]-xbounds[1])/100:xbounds[2]) + xbounds = [A_true - 1.0, A_true + 1.0] + xs = collect(xbounds[1]:((xbounds[2] - xbounds[1]) / 100):xbounds[2]) elseif idx == 3 param = "ω" - xs = collect(ω_true-0.2:0.01:ω_true+0.2) - xbounds = [xs[1], xs[end]] + xs = collect((ω_true - 0.2):0.01:(ω_true + 0.2)) + xbounds = [xs[1], xs[end]] else throw("not implemented") end - if log_normal==true + if log_normal == true xs = log.(xs) - xbounds = log.(xbounds) + xbounds = log.(xbounds) end label = "true " * param posterior_distributions = get_distribution(posterior) - histogram(posterior_distributions[param_names[idx]][1,:], bins=100, normed=true, fill=:slategray, - lab="posterior") + histogram( + posterior_distributions[param_names[idx]][1, :], + bins = 100, + normed = true, + fill = :slategray, + lab = "posterior", + ) prior_plot = get_distribution(mcmc.prior) # This requires StatsPlots - plot!(xs, prior_plot[param_names[idx]], w=2.6, color=:blue, lab="prior") + plot!(xs, prior_plot[param_names[idx]], w = 2.6, color = :blue, lab = "prior") #plot!(xs, mcmc.prior[idx].dist, w=2.6, color=:blue, lab="prior") - plot!([true_values[idx]], seriestype="vline", w=2.6, lab=label) - plot!(xlims=xbounds) + plot!([true_values[idx]], seriestype = "vline", w = 2.6, lab = label) + plot!(xlims = xbounds) title!(param) - + figpath = joinpath(figure_save_directory, "posterior_$(param)_T_$(T)_w_$(ω_true).png") savefig(figpath) linkfig(figpath) end # Save data -@save data_save_directory*"input_output_pairs.jld2" input_output_pairs +@save data_save_directory * "input_output_pairs.jld2" input_output_pairs diff --git a/examples/ci/linkfig.jl b/examples/ci/linkfig.jl index 1e3c5f444..3b08b44cb 100644 --- a/examples/ci/linkfig.jl +++ b/examples/ci/linkfig.jl @@ -1,9 +1,9 @@ function linkfig(figpath) # link figure in logs if we are running on CI if get(ENV, "BUILDKITE", "") == "true" - artifact_url = "artifact://" * join(split(figpath, '/')[end-3:end], '/') + artifact_url = "artifact://" * join(split(figpath, '/')[(end - 3):end], '/') alt = split(splitdir(figpath)[2], '.')[1] @info "Linking Figure: $artifact_url" print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n") end -end \ No newline at end of file +end diff --git a/examples/deprecated/Cloudy/Cloudy_example.jl b/examples/deprecated/Cloudy/Cloudy_example.jl index f4223540c..47e01d763 100644 --- a/examples/deprecated/Cloudy/Cloudy_example.jl +++ b/examples/deprecated/Cloudy/Cloudy_example.jl @@ -92,7 +92,7 @@ dist_true = PDistributions.Gamma(N0_true, θ_true, k_true) ### # Define constraints -lbound_N0 = 0.4 * N0_true +lbound_N0 = 0.4 * N0_true lbound_θ = 1.0e-1 lbound_k = 1.0e-4 c1 = bounded_below(lbound_N0) @@ -125,12 +125,12 @@ n_moments = length(moments) ### # Collision-coalescence kernel to be used in Cloudy -coalescence_coeff = 1/3.14/4/100 +coalescence_coeff = 1 / 3.14 / 4 / 100 kernel_func = x -> coalescence_coeff kernel = Cloudy.KernelTensors.CoalescenceTensor(kernel_func, 0, 100.0) # Time period over which to run Cloudy -tspan = (0., 1.0) +tspan = (0.0, 1.0) ### @@ -140,10 +140,15 @@ tspan = (0., 1.0) ### g_settings_true = GModel.GSettings(kernel, dist_true, moments, tspan) -gt = GModel.run_G(params_true, g_settings_true, PDistributions.update_params, - PDistributions.moment, Cloudy.Sources.get_int_coalescence) +gt = GModel.run_G( + params_true, + g_settings_true, + PDistributions.update_params, + PDistributions.moment, + Cloudy.Sources.get_int_coalescence, +) n_samples = 100 -yt = zeros(length(gt),n_samples) +yt = zeros(length(gt), n_samples) # In a perfect model setting, the "observational noise" represent the internal # model variability. Since Cloudy is a purely deterministic model, there is no # straightforward way of coming up with a covariance structure for this internal @@ -167,9 +172,14 @@ truth_sample = truth.mean N_ens = 50 # number of ensemble members N_iter = 8 # number of EKI iterations # initial parameters: N_params x N_ens -initial_params = EnsembleKalmanProcesses.construct_initial_ensemble(priors, N_ens; rng_seed=6) -ekiobj = EnsembleKalmanProcesses.EnsembleKalmanProcess(initial_params, truth_sample, truth.obs_noise_cov, - Inversion(), Δt=0.1) +initial_params = EnsembleKalmanProcesses.construct_initial_ensemble(priors, N_ens; rng_seed = 6) +ekiobj = EnsembleKalmanProcesses.EnsembleKalmanProcess( + initial_params, + truth_sample, + truth.obs_noise_cov, + Inversion(), + Δt = 0.1, +) # Initialize a ParticleDistribution with dummy parameters. The parameters @@ -180,24 +190,24 @@ g_settings = GModel.GSettings(kernel, dist_type, moments, tspan) # EKI iterations for i in 1:N_iter - - params_i = mapslices(x -> transform_unconstrained_to_constrained(priors, x), - get_u_final(ekiobj); dims=1) - g_ens = GModel.run_G_ensemble(params_i, g_settings, - PDistributions.update_params, - PDistributions.moment, - Cloudy.Sources.get_int_coalescence) + params_i = mapslices(x -> transform_unconstrained_to_constrained(priors, x), get_u_final(ekiobj); dims = 1) + g_ens = GModel.run_G_ensemble( + params_i, + g_settings, + PDistributions.update_params, + PDistributions.moment, + Cloudy.Sources.get_int_coalescence, + ) EnsembleKalmanProcesses.update_ensemble!(ekiobj, g_ens) end # EKI results: Has the ensemble collapsed toward the truth? -transformed_params_true = transform_constrained_to_unconstrained(priors, - params_true) +transformed_params_true = transform_constrained_to_unconstrained(priors, params_true) println("True parameters (transformed): ") println(transformed_params_true) println("\nEKI results:") -println(mean(get_u_final(ekiobj), dims=2)) +println(mean(get_u_final(ekiobj), dims = 2)) ### @@ -207,25 +217,20 @@ println(mean(get_u_final(ekiobj), dims=2)) gppackage = Emulators.GPJL() pred_type = Emulators.YType() gauss_proc = GaussianProcess( - gppackage; - kernel=nothing, # use default squared exponential kernel - prediction_type=pred_type, - noise_learn=false + gppackage; + kernel = nothing, # use default squared exponential kernel + prediction_type = pred_type, + noise_learn = false, ) # Get training points input_output_pairs = Utilities.get_training_points(ekiobj, N_iter) -emulator = Emulator( - gauss_proc, - input_output_pairs, - obs_noise_cov=Γy, - normalize_inputs=true -) +emulator = Emulator(gauss_proc, input_output_pairs, obs_noise_cov = Γy, normalize_inputs = true) optimize_hyperparameters!(emulator) # Check how well the Gaussian Process regression predicts on the # true parameters -y_mean, y_var = Emulators.predict(emulator, transformed_params_true; transform_to_real=true) +y_mean, y_var = Emulators.predict(emulator, reshape(transformed_params_true, :, 1); transform_to_real = true) println("GP prediction on true parameters: ") println(vec(y_mean)) println("true data: ") @@ -237,7 +242,7 @@ println(truth.mean) ### # initial values -u0 = vec(mean(get_inputs(input_output_pairs), dims=2)) +u0 = vec(mean(get_inputs(input_output_pairs), dims = 2)) println("initial parameters: ", u0) # MCMC settings @@ -248,16 +253,14 @@ burnin = 0 step = 0.1 # first guess max_iter = 2000 # number of steps before checking acc/rej rate for step size determination yt_sample = truth_sample -mcmc_test = MarkovChainMonteCarlo.MCMC(yt_sample, Γy, priors, step, u0, max_iter, - mcmc_alg, burnin, svdflag=true) -new_step = MarkovChainMonteCarlo.find_mcmc_step!(mcmc_test, emulator, max_iter=max_iter) +mcmc_test = MarkovChainMonteCarlo.MCMC(yt_sample, Γy, priors, step, u0, max_iter, mcmc_alg, burnin, svdflag = true) +new_step = MarkovChainMonteCarlo.find_mcmc_step!(mcmc_test, emulator, max_iter = max_iter) # Now begin the actual MCMC println("Begin MCMC - with step size ", new_step) burnin = 1000 max_iter = 100000 -mcmc = MarkovChainMonteCarlo.MCMC(yt_sample, Γy, priors, new_step, u0, max_iter, mcmc_alg, - burnin, svdflag=true) +mcmc = MarkovChainMonteCarlo.MCMC(yt_sample, Γy, priors, new_step, u0, max_iter, mcmc_alg, burnin, svdflag = true) MarkovChainMonteCarlo.sample_posterior!(mcmc, emulator, max_iter) posterior = MarkovChainMonteCarlo.get_posterior(mcmc) @@ -273,29 +276,35 @@ println(post_cov) # (in the transformed/unconstrained space) n_params = length(get_name(posterior)) -gr(size=(800,600)) - +gr(size = (800, 600)) + for idx in 1:n_params if idx == 1 - xs = collect(range(5.15, stop=5.25, length=1000)) + xs = collect(range(5.15, stop = 5.25, length = 1000)) elseif idx == 2 - xs = collect(range(0.0, stop=0.5, length=1000)) + xs = collect(range(0.0, stop = 0.5, length = 1000)) elseif idx == 3 - xs = collect(range(-3.0, stop=-2.0, length=1000)) + xs = collect(range(-3.0, stop = -2.0, length = 1000)) else throw("not implemented") end label = "true " * param_names[idx] - posterior_samples = dropdims(get_distribution(posterior)[param_names[idx]], - dims=1) - histogram(posterior_samples, bins=100, normed=true, fill=:slategray, - thickness_scaling=2.0, lab="posterior", legend=:outertopright) + posterior_samples = dropdims(get_distribution(posterior)[param_names[idx]], dims = 1) + histogram( + posterior_samples, + bins = 100, + normed = true, + fill = :slategray, + thickness_scaling = 2.0, + lab = "posterior", + legend = :outertopright, + ) prior_dist = get_distribution(mcmc.prior)[param_names[idx]] - plot!(xs, prior_dist, w=2.6, color=:blue, lab="prior") - plot!([transformed_params_true[idx]], seriestype="vline", w=2.6, lab=label) + plot!(xs, prior_dist, w = 2.6, color = :blue, lab = "prior") + plot!([transformed_params_true[idx]], seriestype = "vline", w = 2.6, lab = label) title!(param_names[idx]) figpath = joinpath(output_directory, "posterior_" * param_names[idx] * ".png") StatsPlots.savefig(figpath) linkfig(figpath) -end \ No newline at end of file +end diff --git a/examples/deprecated/Cloudy/GModel.jl b/examples/deprecated/Cloudy/GModel.jl index d3ff82fee..b8b071303 100644 --- a/examples/deprecated/Cloudy/GModel.jl +++ b/examples/deprecated/Cloudy/GModel.jl @@ -22,53 +22,49 @@ export run_G_ensemble # So, the interface would ideally be something like: # run_G_ensemble(params, G) where G is user-defined """ - GSettings{FT<:AbstractFloat, KT, D} +$(DocStringExtensions.TYPEDEF) -Structure to hold all information to run the forward model G +Structure to hold all information to run the forward model *G*. # Fields -$(DocStringExtensions.FIELDS) +$(DocStringExtensions.TYPEDFIELDS) """ -struct GSettings{FT<:AbstractFloat, KT, D} - "a kernel tensor" +struct GSettings{FT <: AbstractFloat, KT, D} + "A kernel tensor." kernel::KT - "a model distribution" + "A model distribution." dist::D - "the moments of dist that model should return" + "The moments of `dist` that model should return." moments::Array{FT, 1} - "time period over which to run the model, e.g., `(0, 1)`" + "Time period over which to run the model, e.g., `(0, 1)`." tspan::Tuple{FT, FT} end """ - run_G_ensemble(params::Array{FT, 2}, - settings::GSettings{FT}, - update_params, - moment, - get_src; - rng_seed=42) where {FT<:AbstractFloat} - -Run the forward model G for an array of parameters by iteratively -calling run_G for each of the N_ensemble parameter values. -Return g_ens, an array of size N_ensemble x N_data, where -g_ens[j,:] = G(params[j,:]) - - - `params` - array of size N_ensemble x N_parameters containing the - parameters for which G will be run - - `settings` - a GSetttings struct +$(DocStringExtensions.TYPEDSIGNATURES) +Run the forward model *G* for an array of parameters by iteratively +calling `run_G` for each of the *N\\_ensemble* parameter values. + +- `params` - array of size (*N\\_ensemble* × *N\\_parameters*) containing the parameters for which + G will be run. +- `settings` - a [GSetttings](@ref) struct. + +Returns `g_ens``, an array of size (*N\\_ensemble* × *N\\_data*), where g_ens[j,:] = G(params[j,:]). """ -function run_G_ensemble(params::Array{FT, 2}, - settings::GSettings{FT}, - update_params, - moment, - get_src; - rng_seed=42) where {FT<:AbstractFloat} +function run_G_ensemble( + params::Array{FT, 2}, + settings::GSettings{FT}, + update_params, + moment, + get_src; + rng_seed = 42, +) where {FT <: AbstractFloat} N_ens = size(params, 2) # params is N_ens x N_params n_moments = length(settings.moments) - g_ens = zeros(n_moments,N_ens) + g_ens = zeros(n_moments, N_ens) Random.seed!(rng_seed) for i in 1:N_ens @@ -81,23 +77,14 @@ end """ - run_G(u::Array{FT, 1}, - settings::GSettings{FT}, - update_params, - moment, - get_src) where {FT<:AbstractFloat} - -Return g_u = G(u), a vector of length N_data. +$(DocStringExtensions.TYPEDSIGNATURES) - - `u` - parameter vector of length N_parameters - - `settings` - a GSettings struct +- `u` - parameter vector of length *N\\_parameters*. +- `settings` - a [GSetttings](@ref) struct. +Returns `g_u` = *G(u)*, a vector of length *N\\_data*. """ -function run_G(u::Array{FT, 1}, - settings::GSettings{FT}, - update_params, - moment, - get_src) where {FT<:AbstractFloat} +function run_G(u::Array{FT, 1}, settings::GSettings{FT}, update_params, moment, get_src) where {FT <: AbstractFloat} # generate the initial distribution dist = update_params(settings.dist, u) @@ -117,7 +104,7 @@ function run_G(u::Array{FT, 1}, rhs(M, p, t) = get_src(M, dist, settings.kernel) prob = ODEProblem(rhs, moments_init, settings.tspan) # Solve the ODE - sol = solve(prob, CVODE_BDF(), alg_hints=[:stiff], reltol=tol, abstol=tol) + sol = solve(prob, CVODE_BDF(), alg_hints = [:stiff], reltol = tol, abstol = tol) # Return moments at last time step moments_final = vcat(sol.u'...)[end, :] diff --git a/src/CalibrateEmulateSample.jl b/src/CalibrateEmulateSample.jl index 9a233e192..9f227d0ee 100644 --- a/src/CalibrateEmulateSample.jl +++ b/src/CalibrateEmulateSample.jl @@ -1,17 +1,18 @@ +""" +# Imported modules: +$(IMPORTS) + +# Exports: +$(EXPORTS) +""" module CalibrateEmulateSample using Distributions, Statistics, LinearAlgebra, DocStringExtensions # imported modules from EKP. -import EnsembleKalmanProcesses: EnsembleKalmanProcesses, - ParameterDistributions, - Observations, - DataContainers - -export EnsembleKalmanProcesses, - ParameterDistributions, - Observations, - DataContainers +import EnsembleKalmanProcesses: EnsembleKalmanProcesses, ParameterDistributions, Observations, DataContainers + +export EnsembleKalmanProcesses, ParameterDistributions, Observations, DataContainers # No internal deps, heavy external deps #include("GaussianProcessEmulator.jl") diff --git a/src/Emulator.jl b/src/Emulator.jl index 8b6f12873..f7bc3bf4b 100644 --- a/src/Emulator.jl +++ b/src/Emulator.jl @@ -1,6 +1,8 @@ module Emulators using ..DataContainers + +using DocStringExtensions using Statistics using Distributions using LinearAlgebra @@ -12,7 +14,7 @@ export optimize_hyperparameters! export predict """ - MachineLearningTool +$(DocStringExtensions.TYPEDEF) Type to dispatch different emulators: @@ -23,9 +25,15 @@ abstract type MachineLearningTool end function throw_define_mlt() throw(ErrorException("Unknown MachineLearningTool defined, please use a known implementation")) end -function build_models!(mlt,iopairs) throw_define_mlt() end -function optimize_hyperparameters!(mlt) throw_define_mlt() end -function predict(mlt,new_inputs) throw_define_mlt() end +function build_models!(mlt, iopairs) + throw_define_mlt() +end +function optimize_hyperparameters!(mlt) + throw_define_mlt() +end +function predict(mlt, new_inputs) + throw_define_mlt() +end # include the different <: ML models include("GaussianProcess.jl") #for GaussianProcess @@ -37,26 +45,29 @@ include("GaussianProcess.jl") #for GaussianProcess # We will define the different emulator types after the general statements """ - Emulator +$(DocStringExtensions.TYPEDEF) -Structure used to represent a general emulator: +Structure used to represent a general emulator, independently of the algorithm used. + +# Fields +$(DocStringExtensions.TYPEDFIELDS) """ -struct Emulator{FT<:AbstractFloat} - "Machine learning tool, defined as a struct of type MachineLearningTool" +struct Emulator{FT <: AbstractFloat} + "Machine learning tool, defined as a struct of type MachineLearningTool." machine_learning_tool::MachineLearningTool - "normalized, standardized, transformed pairs given the Boolean's normalize_inputs, standardize_outputs, truncate_svd " + "Normalized, standardized, transformed pairs given the Boolean's normalize\\_inputs, standardize\\_outputs, truncate\\_svd." training_pairs::PairedDataContainer{FT} - "mean of input; length input_dim" + "Mean of input; length *input\\_dim*." input_mean::AbstractVector{FT} - "square root of the inverse of the input covariance matrix; input_dim × input_dim" + "Square root of the inverse of the input covariance matrix; size *input\\_dim* × *input\\_dim*." normalize_inputs::Bool - "whether to fit models on normalized outputs outputs / standardize_outputs_factor" + "Whether to fit models on normalized outputs: `outputs / standardize_outputs_factor`." sqrt_inv_input_cov::Union{AbstractMatrix{FT}, UniformScaling{FT}, Nothing} - " if normalizing: whether to fit models on normalized inputs ((inputs - input_mean) * sqrt_inv_input_cov)" + "If normalizing: whether to fit models on normalized inputs (`(inputs - input_mean) * sqrt_inv_input_cov`)." standardize_outputs::Bool - "if standardizing: Standardization factors (characteristic values of the problem)" + "If standardizing: Standardization factors (characteristic values of the problem)." standardize_outputs_factors::Union{AbstractVector{FT}, Nothing} - "the singular value decomposition of obs_noise_cov, such that obs_noise_cov = decomposition.U * Diagonal(decomposition.S) * decomposition.Vt. NB: the svd may be reduced in dimensions" + "The singular value decomposition of *obs\\_noise\\_cov*, such that *obs\\_noise\\_cov* = decomposition.U * Diagonal(decomposition.S) * decomposition.Vt. NB: the SVD may be reduced in dimensions." decomposition::Union{SVD, Nothing} end @@ -68,8 +79,8 @@ function Emulator( normalize_inputs::Bool = true, standardize_outputs::Bool = false, standardize_outputs_factors::Union{AbstractVector{FT}, Nothing} = nothing, - truncate_svd::FT=1.0 -) where {FT<:AbstractFloat} + truncate_svd::FT = 1.0, +) where {FT <: AbstractFloat} # For Consistency checks input_dim, output_dim = size(input_output_pairs, 1) @@ -83,41 +94,35 @@ function Emulator( size(obs_noise_cov) == (output_dim, output_dim) || throw(ArgumentError(err2)) end - + # [1.] Normalize the inputs? - input_mean = vec(mean(get_inputs(input_output_pairs), dims=2)) #column vector + input_mean = vec(mean(get_inputs(input_output_pairs), dims = 2)) #column vector sqrt_inv_input_cov = nothing if normalize_inputs # Normalize (NB the inputs have to be of) size [input_dim × N_samples] to pass to ML tool - sqrt_inv_input_cov = sqrt(inv(Symmetric(cov(get_inputs(input_output_pairs), dims=2)))) - training_inputs = normalize( - get_inputs(input_output_pairs), - input_mean, - sqrt_inv_input_cov - ) + sqrt_inv_input_cov = sqrt(inv(Symmetric(cov(get_inputs(input_output_pairs), dims = 2)))) + training_inputs = normalize(get_inputs(input_output_pairs), input_mean, sqrt_inv_input_cov) else training_inputs = get_inputs(input_output_pairs) end # [2.] Standardize the outputs? if standardize_outputs - training_outputs, obs_noise_cov = standardize( - get_outputs(input_output_pairs), - obs_noise_cov, - standardize_outputs_factors - ) + training_outputs, obs_noise_cov = + standardize(get_outputs(input_output_pairs), obs_noise_cov, standardize_outputs_factors) else - training_outputs = get_outputs(input_output_pairs) + training_outputs = get_outputs(input_output_pairs) end # [3.] Decorrelating the outputs [always performed] #Transform data if obs_noise_cov available # (if obs_noise_cov==nothing, transformed_data is equal to data) - decorrelated_training_outputs, decomposition = svd_transform(training_outputs, obs_noise_cov, truncate_svd=truncate_svd) + decorrelated_training_outputs, decomposition = + svd_transform(training_outputs, obs_noise_cov, truncate_svd = truncate_svd) # write new pairs structure - if truncate_svd<1.0 + if truncate_svd < 1.0 #note this changes the dimension of the outputs training_pairs = PairedDataContainer(training_inputs, decorrelated_training_outputs) input_dim, output_dim = size(training_pairs, 1) @@ -126,8 +131,8 @@ function Emulator( end # [4.] build an emulator - build_models!(machine_learning_tool,training_pairs) - + build_models!(machine_learning_tool, training_pairs) + return Emulator{FT}( machine_learning_tool, training_pairs, @@ -136,39 +141,41 @@ function Emulator( sqrt_inv_input_cov, standardize_outputs, standardize_outputs_factors, - decomposition + decomposition, ) -end +end """ - function optimize_hyperparameters!(emulator::Emulator) +$(DocStringExtensions.TYPEDSIGNATURES) -optimize the hyperparameters in the machine learning tool +Optimizes the hyperparameters in the machine learning tool. """ -function optimize_hyperparameters!(emulator::Emulator{FT}) where {FT<:AbstractFloat} +function optimize_hyperparameters!(emulator::Emulator{FT}) where {FT <: AbstractFloat} optimize_hyperparameters!(emulator.machine_learning_tool) end """ - function predict(emulator::Emulator, new_inputs; transform_to_real=false) +$(DocStringExtensions.TYPEDSIGNATURES) -makes a prediction using the emulator on new inputs (each new inputs given as data columns), default is to predict in the decorrelated space +Makes a prediction using the emulator on new inputs (each new inputs given as data columns). +Default is to predict in the decorrelated space. """ function predict( - emulator::Emulator{FT}, - new_inputs::AbstractMatrix{FT}; - transform_to_real=false -) where {FT<:AbstractFloat} + emulator::Emulator{FT}, + new_inputs::AbstractMatrix{FT}; + transform_to_real = false, +) where {FT <: AbstractFloat} # Check if the size of new_inputs is consistent with the GP model's input # dimension. input_dim, output_dim = size(emulator.training_pairs, 1) N_samples = size(new_inputs, 2) - size(new_inputs, 1) == input_dim || throw(ArgumentError("GP object and input observations do not have consistent dimensions")) + size(new_inputs, 1) == input_dim || + throw(ArgumentError("GP object and input observations do not have consistent dimensions")) # [1.] normalize - normalized_new_inputs = normalize(emulator,new_inputs) + normalized_new_inputs = normalize(emulator, new_inputs) # [2.] predict. Note: ds = decorrelated, standard ds_outputs, ds_output_var = predict(emulator.machine_learning_tool, normalized_new_inputs) @@ -180,14 +187,13 @@ function predict( Try setting transform_to_real=false""")) elseif transform_to_real && emulator.decomposition !== nothing #transform back to real coords - cov becomes dense - s_outputs, s_output_cov = svd_reverse_transform_mean_cov( - ds_outputs, ds_output_var, emulator.decomposition) + s_outputs, s_output_cov = svd_reverse_transform_mean_cov(ds_outputs, ds_output_var, emulator.decomposition) if output_dim == 1 s_output_cov = [s_output_cov[i][1] for i in 1:N_samples] end # [4.] unstandardize return reverse_standardize(emulator, s_outputs, s_output_cov) - + else # remain in decorrelated, standardized coordinates (cov remains diagonal) # Convert to vector of matrices to match the format @@ -198,7 +204,7 @@ function predict( end return ds_outputs, ds_output_diagvar - + end end @@ -208,11 +214,11 @@ end # Normalization, Standardization, and Decorrelation """ - function normalize(emulator::Emulator, inputs) +$(DocStringExtensions.TYPEDSIGNATURES) -normalize the input data, with a normalizing function +Normalize the input data, with a normalizing function. """ -function normalize(emulator::Emulator{FT}, inputs::AbstractVecOrMat{FT}) where {FT<:AbstractFloat} +function normalize(emulator::Emulator{FT}, inputs::AbstractVecOrMat{FT}) where {FT <: AbstractFloat} if emulator.normalize_inputs return normalize(inputs, emulator.input_mean, emulator.sqrt_inv_input_cov) else @@ -221,29 +227,29 @@ function normalize(emulator::Emulator{FT}, inputs::AbstractVecOrMat{FT}) where { end """ - function normalize(inputs, input_mean, sqrt_inv_input_cov) +$(DocStringExtensions.TYPEDSIGNATURES) -normalize with the empirical Gaussian distribution of points +Normalize with the empirical Gaussian distribution of points. """ function normalize( - inputs::AbstractVecOrMat{FT}, - input_mean::AbstractVector{FT}, - sqrt_inv_input_cov::Union{AbstractMatrix{FT}, UniformScaling{FT}} -) where {FT<:AbstractFloat} + inputs::AbstractVecOrMat{FT}, + input_mean::AbstractVector{FT}, + sqrt_inv_input_cov::Union{AbstractMatrix{FT}, UniformScaling{FT}}, +) where {FT <: AbstractFloat} training_inputs = sqrt_inv_input_cov * (inputs .- input_mean) - return training_inputs + return training_inputs end """ - function standardize(outputs, output_cov, factors) +$(DocStringExtensions.TYPEDSIGNATURES) -standardize with a vector of factors (size equal to output dimension) +Standardize with a vector of factors (size equal to output dimension). """ function standardize( - outputs::AbstractVecOrMat{FT}, - output_covs::Vector{<:Union{AbstractMatrix{FT}, UniformScaling{FT}}}, - factors::AbstractVector{FT} -) where {FT<:AbstractFloat} + outputs::AbstractVecOrMat{FT}, + output_covs::Vector{<:Union{AbstractMatrix{FT}, UniformScaling{FT}}}, + factors::AbstractVector{FT}, +) where {FT <: AbstractFloat} # Case where `output_covs` is a Vector of covariance matrices n = length(factors) outputs = outputs ./ factors @@ -254,8 +260,8 @@ function standardize( elseif isa(output_covs[i], UniformScaling) # assert all elements of factors are same, otherwise can't cast back to # UniformScaling - @assert all(y -> y == first(factors), factors) - output_covs[i] = output_covs[i] / cov_factors[1,1] + @assert all(y -> y == first(factors), factors) + output_covs[i] = output_covs[i] / cov_factors[1, 1] else # Diagonal, TriDiagonal etc. case # need to do ./ as dense matrix and cast back to original type @@ -269,27 +275,27 @@ function standardize( end function standardize( - outputs::AbstractVecOrMat{FT}, - output_cov::Union{AbstractMatrix{FT}, UniformScaling{FT}}, - factors::AbstractVector{FT} -) where {FT<:AbstractFloat} + outputs::AbstractVecOrMat{FT}, + output_cov::Union{AbstractMatrix{FT}, UniformScaling{FT}}, + factors::AbstractVector{FT}, +) where {FT <: AbstractFloat} # Case where `output_cov` is a single covariance matrix stdized_out, stdized_covs = standardize(outputs, [output_cov], factors) return stdized_out, stdized_covs[1] end """ - function reverse_standardize(emulator::Emulator, outputs, output_covs) +$(DocStringExtensions.TYPEDSIGNATURES) -reverse a previous standardization with the stored vector of factors (size equal to output +Reverse a previous standardization with the stored vector of factors (size equal to output dimension). `output_cov` is a Vector of covariance matrices, such as is returned by [`svd_reverse_transform_mean_cov`](@ref). """ function reverse_standardize( - emulator::Emulator{FT}, - outputs::AbstractVecOrMat{FT}, - output_covs::Union{AbstractMatrix{FT}, Vector{<:AbstractMatrix{FT}}} -) where {FT<:AbstractFloat} + emulator::Emulator{FT}, + outputs::AbstractVecOrMat{FT}, + output_covs::Union{AbstractMatrix{FT}, Vector{<:AbstractMatrix{FT}}}, +) where {FT <: AbstractFloat} if emulator.standardize_outputs return standardize(outputs, output_covs, 1.0 ./ emulator.standardize_outputs_factors) else @@ -300,10 +306,10 @@ end """ - svd_transform(data, obs_noise_cov, truncate_svd) where {FT} +$(DocStringExtensions.TYPEDSIGNATURES) Apply a singular value decomposition (SVD) to the data - - `data` - GP training data/targets; output_dim × N_samples + - `data` - GP training data/targets; size *output\\_dim* × *N\\_samples* - `obs_noise_cov` - covariance of observational noise - `truncate_svd` - Project onto this fraction of the largest principal components. Defaults to 1.0 (no truncation). @@ -311,23 +317,23 @@ Apply a singular value decomposition (SVD) to the data Returns the transformed data and the decomposition, which is a matrix factorization of type LinearAlgebra.SVD. -Note: If F::SVD is the factorization object, U, S, V and Vt can be obtained via +Note: If `F::SVD` is the factorization object, U, S, V and Vt can be obtained via F.U, F.S, F.V and F.Vt, such that A = U * Diagonal(S) * Vt. The singular values in S are sorted in descending order. """ function svd_transform( - data::AbstractMatrix{FT}, - obs_noise_cov::Union{AbstractMatrix{FT}, Nothing}; - truncate_svd::FT=1.0 -) where {FT<:AbstractFloat} + data::AbstractMatrix{FT}, + obs_noise_cov::Union{AbstractMatrix{FT}, Nothing}; + truncate_svd::FT = 1.0, +) where {FT <: AbstractFloat} if obs_noise_cov === nothing # no-op case return data, nothing end # actually have a matrix to take the SVD of decomp = svd(obs_noise_cov) - sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomp.S)) - if truncate_svd < 1.0 + sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomp.S)) + if truncate_svd < 1.0 # Truncate the SVD as a form of regularization # Find cutoff S_cumsum = cumsum(decomp.S) / sum(decomp.S) @@ -336,70 +342,69 @@ function svd_transform( n = size(data)[1] println("SVD truncated at k: ", k, "/", n) # Apply truncated SVD - transformed_data = sqrt_singular_values_inv[1:k, 1:k] * decomp.Vt[1:k, :] * data + transformed_data = sqrt_singular_values_inv[1:k, 1:k] * decomp.Vt[1:k, :] * data decomposition = SVD(decomp.U[:, 1:k], decomp.S[1:k], decomp.Vt[1:k, :]) - else + else transformed_data = sqrt_singular_values_inv * decomp.Vt * data - decomposition = svd(obs_noise_cov) + decomposition = svd(obs_noise_cov) end return transformed_data, decomposition end function svd_transform( - data::AbstractVector{FT}, - obs_noise_cov::Union{AbstractMatrix{FT}, Nothing}; - truncate_svd::FT = 1.0 -) where {FT<:AbstractFloat} + data::AbstractVector{FT}, + obs_noise_cov::Union{AbstractMatrix{FT}, Nothing}; + truncate_svd::FT = 1.0, +) where {FT <: AbstractFloat} # method for 1D data - transformed_data, decomposition = svd_transform( - reshape(data, :, 1), obs_noise_cov; truncate_svd=truncate_svd - ) + transformed_data, decomposition = svd_transform(reshape(data, :, 1), obs_noise_cov; truncate_svd = truncate_svd) return vec(transformed_data), decomposition end function svd_transform( - data::AbstractVecOrMat{FT}, obs_noise_cov::UniformScaling{FT}; - truncate_svd::FT = 1.0 -) where {FT<:AbstractFloat} + data::AbstractVecOrMat{FT}, + obs_noise_cov::UniformScaling{FT}; + truncate_svd::FT = 1.0, +) where {FT <: AbstractFloat} # method for UniformScaling - return svd_transform(data, Diagonal(obs_noise_cov, size(data, 1)); truncate_svd=truncate_svd) + return svd_transform(data, Diagonal(obs_noise_cov, size(data, 1)); truncate_svd = truncate_svd) end """ -svd_reverse_transform_mean_cov(μ, σ2, decomposition::SVD) where {FT} +$(DocStringExtensions.TYPEDSIGNATURES) Transform the mean and covariance back to the original (correlated) coordinate system - - `μ` - predicted mean; output_dim × N_predicted_points - - `σ2` - predicted variance; output_dim × N_predicted_points - - `decomposition` - SVD decomposition of obs_noise_cov. + - `μ` - predicted mean; size *output\\_dim* × *N\\_predicted\\_points*. + - `σ2` - predicted variance; size *output\\_dim* × *N\\_predicted\\_points*. + - `decomposition` - SVD decomposition of *obs\\_noise\\_cov*. -Returns the transformed mean (output_dim × N_predicted_points) and variance. +Returns the transformed mean (size *output\\_dim* × *N\\_predicted\\_points*) and variance. Note that transforming the variance back to the original coordinate system results in non-zero off-diagonal elements, so instead of just returning the elements on the main diagonal (i.e., the variances), we return the full -covariance at each point, as a vector of length N_predicted_points, where -each element is a matrix of size output_dim × output_dim +covariance at each point, as a vector of length *N\\_predicted\\_points*, where +each element is a matrix of size *output\\_dim* × *output\\_dim*. """ function svd_reverse_transform_mean_cov( μ::AbstractMatrix{FT}, - σ2::AbstractMatrix{FT}, - decomposition::SVD -) where {FT<:AbstractFloat} + σ2::AbstractMatrix{FT}, + decomposition::SVD, +) where {FT <: AbstractFloat} @assert size(μ) == size(σ2) - + output_dim, N_predicted_points = size(σ2) # We created meanvGP = D_inv * Vt * mean_v so meanv = V * D * meanvGP sqrt_singular_values = Diagonal(sqrt.(decomposition.S)) transformed_μ = decomposition.V * sqrt_singular_values * μ - + transformed_σ2 = [zeros(output_dim, output_dim) for i in 1:N_predicted_points] # Back transformation - + for j in 1:N_predicted_points - σ2_j = decomposition.V * sqrt_singular_values * Diagonal(σ2[:,j]) * sqrt_singular_values * decomposition.Vt + σ2_j = decomposition.V * sqrt_singular_values * Diagonal(σ2[:, j]) * sqrt_singular_values * decomposition.Vt transformed_σ2[j] = σ2_j end - + return transformed_μ, transformed_σ2 end diff --git a/src/GaussianProcess.jl b/src/GaussianProcess.jl index afb4b2880..9662385d4 100644 --- a/src/GaussianProcess.jl +++ b/src/GaussianProcess.jl @@ -2,7 +2,8 @@ using GaussianProcesses using EnsembleKalmanProcesses.DataContainers - + +using DocStringExtensions using PyCall using ScikitLearn const pykernels = PyNULL() @@ -19,95 +20,88 @@ export GPJL, SKLJL export YType, FType """ - GaussianProcessesPackage +$(DocStringExtensions.TYPEDEF) -Type to dispatch which GP package to use +Type to dispatch which GP package to use: - `GPJL` for GaussianProcesses.jl, - - `SKLJL` for the ScikitLearn GaussianProcessRegressor + - `SKLJL` for the ScikitLearn GaussianProcessRegressor. """ abstract type GaussianProcessesPackage end struct GPJL <: GaussianProcessesPackage end struct SKLJL <: GaussianProcessesPackage end """ - PredictionType +$(DocStringExtensions.TYPEDEF) -Predict type for `GPJL` in GaussianProcesses.jl +Predict type for `GPJL` in GaussianProcesses.jl: - `YType` - - `FType` latent function + - `FType` latent function. """ abstract type PredictionType end struct YType <: PredictionType end struct FType <: PredictionType end """ - GaussianProcess{FT<:AbstractFloat} <: MachineLearningTool +$(DocStringExtensions.TYPEDEF) - Structure holding training input and the fitted Gaussian Process Regression +Structure holding training input and the fitted Gaussian process regression models. # Fields -$(DocStringExtensions.FIELDS) +$(DocStringExtensions.TYPEDFIELDS) """ struct GaussianProcess{GPPackage} <: MachineLearningTool - "the Gaussian Process (GP) Regression model(s) that are fitted to the given input-data pairs" + "The Gaussian Process (GP) Regression model(s) that are fitted to the given input-data pairs." models::Vector{Union{<:GaussianProcesses.GPE, <:PyObject, Nothing}} - "Kernel object" - kernel::Union{<:Kernel, <:PyObject, Nothing} - "learn the noise with the White Noise kernel explicitly?" + "Kernel object." + kernel::Union{<:Kernel, <:PyObject, Nothing} + "Learn the noise with the White Noise kernel explicitly?" noise_learn::Bool - "prediction type (`y` to predict the data, `f` to predict the latent function)" + "Prediction type (`y` to predict the data, `f` to predict the latent function)." prediction_type::PredictionType end """ - function GaussianProcess( - package; - kernel::Union{K, KPy, Nothing}=nothing, - noise_learn=true, - prediction_type::PredictionType=YType(), - ) where {K<:Kernel, KPy<:PyObject} - - - `package` - GaussianProcessPackage object - - `kernel` - GaussianProcesses kernel object. default is a Squared Exponential kernel. - - `noise_learn` - boolean to additionally learn white noise in decorrelated space. default is true. - - `prediction_type` - PredictionType object. default predicts data, not latent function (FType()). +$(DocStringExtensions.TYPEDSIGNATURES) + + - `package` - GaussianProcessPackage object. + - `kernel` - GaussianProcesses kernel object. Default is a Squared Exponential kernel. + - `noise_learn` - Boolean to additionally learn white noise in decorrelated space. Default is true. + - `prediction_type` - PredictionType object. Default predicts data, not latent function (FType()). """ function GaussianProcess( package::GPPkg; kernel::Union{K, KPy, Nothing} = nothing, noise_learn = true, prediction_type::PredictionType = YType(), -) where {GPPkg<:GaussianProcessesPackage, K<:Kernel, KPy<:PyObject} +) where {GPPkg <: GaussianProcessesPackage, K <: Kernel, KPy <: PyObject} # Initialize vector for GP models models = Vector{Union{<:GaussianProcesses.GPE, <:PyObject, Nothing}}(undef, 0) - + return GaussianProcess{typeof(package)}(models, kernel, noise_learn, prediction_type) end # First we create the GPJL implementation """ - function build_models!( - gp::GaussianProcess{package}, - input_output_pairs) +$(DocStringExtensions.TYPEDSIGNATURES) -method to build gaussian process models based on the package +Method to build Gaussian process models based on the package. """ function build_models!( gp::GaussianProcess{GPJL}, - input_output_pairs::PairedDataContainer{FT} -) where {FT<:AbstractFloat} + input_output_pairs::PairedDataContainer{FT}, +) where {FT <: AbstractFloat} # get inputs and outputs input_values = get_inputs(input_output_pairs) output_values = get_outputs(input_output_pairs) # Number of models (We are fitting one model per output dimension, as data is decorrelated) models = gp.models - N_models = size(output_values,1); #size(transformed_data)[1] + N_models = size(output_values, 1) #size(transformed_data)[1] # Use a default kernel unless a kernel was supplied to GaussianProcess @@ -132,57 +126,53 @@ function build_models!( white = Noise(white_logstd) kern = kern + white println("Learning additive white noise") - + # Make the regularization small. We actually learn # total_noise = white_logstd + logstd_regularization_noise magic_number = 1e-3 # magic_number << 1 - logstd_regularization_noise = log(sqrt(magic_number)) + logstd_regularization_noise = log(sqrt(magic_number)) else # When not learning noise, our SVD transformation implies the # observational noise is the identity. - logstd_regularization_noise = log(sqrt(1.0)) + logstd_regularization_noise = log(sqrt(1.0)) end - + for i in 1:N_models # Make a copy of the kernel (because it gets altered in every # iteration) kernel_i = deepcopy(kern) println("kernel in GaussianProcess:") println(kernel_i) - data_i = output_values[i,:] + data_i = output_values[i, :] # GaussianProcesses.GPE() arguments: # input_values: (input_dim × N_samples) # GPdata_i: (N_samples,) - + # Zero mean function kmean = MeanZero() # Instantiate GP model - m = GaussianProcesses.GPE(input_values, - output_values[i,:], - kmean, - kernel_i, - logstd_regularization_noise) - + m = GaussianProcesses.GPE(input_values, output_values[i, :], kmean, kernel_i, logstd_regularization_noise) + println("created GP: ", i) push!(models, m) - + end end """ - function optimize_hyperparameters!(gp::GaussianProcess{package}) +$(DocStringExtensions.TYPEDSIGNATURES) -optimize Gaussian Process hyperparameters using in-build package method +Optimize Gaussian process hyperparameters using in-build package method. """ function optimize_hyperparameters!(gp::GaussianProcess{GPJL}) N_models = length(gp.models) - for i = 1:N_models + for i in 1:N_models # always regress with noise_learn=false; if gp was created with noise_learn=true # we've already explicitly added noise to the kernel - optimize!(gp.models[i], noise=false) + optimize!(gp.models[i], noise = false) println("optimized hyperparameters of GP: ", i) println(gp.models[i].kernel) end @@ -190,58 +180,56 @@ end # subroutine with common predict() logic function _predict( - gp::GaussianProcess, - new_inputs::AbstractMatrix{FT}, - predict_method::Function -) where {FT<:AbstractFloat} + gp::GaussianProcess, + new_inputs::AbstractMatrix{FT}, + predict_method::Function, +) where {FT <: AbstractFloat} M = length(gp.models) N_samples = size(new_inputs, 2) # Predicts columns of inputs: input_dim × N_samples μ = zeros(M, N_samples) σ2 = zeros(M, N_samples) for i in 1:M - μ[i,:], σ2[i,:] = predict_method(gp.models[i], new_inputs) + μ[i, :], σ2[i, :] = predict_method(gp.models[i], new_inputs) end return μ, σ2 end -predict(gp::GaussianProcess{GPJL}, new_inputs::AbstractMatrix{FT}, ::YType) where {FT<:AbstractFloat} = +predict(gp::GaussianProcess{GPJL}, new_inputs::AbstractMatrix{FT}, ::YType) where {FT <: AbstractFloat} = _predict(gp, new_inputs, GaussianProcesses.predict_y) -predict(gp::GaussianProcess{GPJL}, new_inputs::AbstractMatrix{FT}, ::FType) where {FT<:AbstractFloat} = +predict(gp::GaussianProcess{GPJL}, new_inputs::AbstractMatrix{FT}, ::FType) where {FT <: AbstractFloat} = _predict(gp, new_inputs, GaussianProcesses.predict_f) """ - function predict(gp::GaussianProcess{package}, new_inputs) +$(DocStringExtensions.TYPEDSIGNATURES) -predict means and covariances in decorrelated output space using gaussian process models +Predict means and covariances in decorrelated output space using Gaussian process models. """ -predict(gp::GaussianProcess{GPJL}, new_inputs::AbstractMatrix{FT}) where {FT<:AbstractFloat} = +predict(gp::GaussianProcess{GPJL}, new_inputs::AbstractMatrix{FT}) where {FT <: AbstractFloat} = predict(gp, new_inputs, gp.prediction_type) #now we build the SKLJL implementation function build_models!( gp::GaussianProcess{SKLJL}, - input_output_pairs::PairedDataContainer{FT} -) where {FT<:AbstractFloat} - # get inputs and outputs - input_values = permutedims(get_inputs(input_output_pairs), (2,1)) + input_output_pairs::PairedDataContainer{FT}, +) where {FT <: AbstractFloat} + # get inputs and outputs + input_values = permutedims(get_inputs(input_output_pairs), (2, 1)) output_values = get_outputs(input_output_pairs) # Number of models (We are fitting one model per output dimension, as data is decorrelated) models = gp.models - N_models = size(output_values,1); #size(transformed_data)[1] + N_models = size(output_values, 1) #size(transformed_data)[1] if gp.kernel === nothing println("Using default squared exponential kernel, learning length scale and variance parameters") # Create default squared exponential kernel const_value = 1.0 - var_kern = pykernels.ConstantKernel(constant_value=const_value, - constant_value_bounds=(1e-5, 1e4)) + var_kern = pykernels.ConstantKernel(constant_value = const_value, constant_value_bounds = (1e-5, 1e4)) rbf_len = ones(size(input_values, 2)) - rbf = pykernels.RBF(length_scale=rbf_len, - length_scale_bounds=(1e-5, 1e4)) + rbf = pykernels.RBF(length_scale = rbf_len, length_scale_bounds = (1e-5, 1e4)) kern = var_kern * rbf println("Using default squared exponential kernel:", kern) else @@ -252,8 +240,7 @@ function build_models!( if gp.noise_learn # Add white noise to kernel white_noise_level = 1.0 - white = pykernels.WhiteKernel(noise_level=white_noise_level, - noise_level_bounds=(1e-05, 10.0)) + white = pykernels.WhiteKernel(noise_level = white_noise_level, noise_level_bounds = (1e-05, 10.0)) kern = kern + white println("Learning additive white noise") @@ -266,19 +253,17 @@ function build_models!( # observational noise is the identity. regularization_noise = 1.0 end - + for i in 1:N_models kernel_i = deepcopy(kern) data_i = output_values[i, :] - m = pyGP.GaussianProcessRegressor(kernel=kernel_i, - n_restarts_optimizer=10, - alpha=regularization_noise) + m = pyGP.GaussianProcessRegressor(kernel = kernel_i, n_restarts_optimizer = 10, alpha = regularization_noise) # ScikitLearn.fit! arguments: # input_values: (N_samples × input_dim) # data_i: (N_samples,) ScikitLearn.fit!(m, input_values, data_i) - if i==1 + if i == 1 println(m.kernel.hyperparameters) print("Completed training of: ") end @@ -293,11 +278,10 @@ function optimize_hyperparameters!(gp::GaussianProcess{SKLJL}) println("SKlearn, already trained. continuing...") end -function _SKJL_predict_function(gp_model::PyObject, new_inputs::AbstractMatrix{FT}) where {FT<:AbstractFloat} +function _SKJL_predict_function(gp_model::PyObject, new_inputs::AbstractMatrix{FT}) where {FT <: AbstractFloat} # SKJL based on rows not columns; need to transpose inputs - μ, σ = gp_model.predict(new_inputs', return_std=true) + μ, σ = gp_model.predict(new_inputs', return_std = true) return μ, (σ .* σ) end -predict(gp::GaussianProcess{SKLJL}, new_inputs::AbstractMatrix{FT}) where {FT<:AbstractFloat} = +predict(gp::GaussianProcess{SKLJL}, new_inputs::AbstractMatrix{FT}) where {FT <: AbstractFloat} = _predict(gp, new_inputs, _SKJL_predict_function) - diff --git a/src/MCMC_sampler_defaults.jl b/src/MCMC_sampler_defaults.jl new file mode 100644 index 000000000..cf6823bcd --- /dev/null +++ b/src/MCMC_sampler_defaults.jl @@ -0,0 +1,251 @@ +# Default methods (and other boilerplate) for the AbstractMCMC interface. +# +# "Default" here means the logic in this file is independent of the specific samplers +# implemented in MCMC_samplers.jl. +# Code here is taken directly from AdvancedMH.jl; we don't simply add that package as a +# dependency because the fraction of code we need is relatively small. + + +# ------------------------------------------------------------------------------------------ +# Associated methods for Proposal objects + +# Random draws +Base.rand(p::Proposal, args...) = rand(Random.GLOBAL_RNG, p, args...) +Base.rand(rng::Random.AbstractRNG, p::Proposal{<:Distribution}) = rand(rng, p.proposal) +function Base.rand(rng::Random.AbstractRNG, p::Proposal{<:AbstractArray}) + return map(x -> rand(rng, x), p.proposal) +end + +# Densities +Distributions.logpdf(p::Proposal{<:Distribution}, v) = logpdf(p.proposal, v) +function Distributions.logpdf(p::Proposal{<:AbstractArray}, v) + # `mapreduce` with multiple iterators requires Julia 1.2 or later + return mapreduce(((pi, vi),) -> logpdf(pi, vi), +, zip(p.proposal, v)) +end + +# Proposals +function propose(rng::Random.AbstractRNG, proposal::Proposal{<:Function}, model::DensityModel) + return propose(rng, proposal(), model) +end + +function propose(rng::Random.AbstractRNG, proposal::Proposal{<:Function}, model::DensityModel, t) + return propose(rng, proposal(t), model) +end + +function q(proposal::Proposal{<:Function}, t, t_cond) + return q(proposal(t_cond), t, t_cond) +end + +# ------------------------------------------------------------------------------------------ +# Extension of propose() method to multiple proposals + +function propose(rng::Random.AbstractRNG, proposals::AbstractArray{<:Proposal}, model::DensityModel) + return map(proposals) do proposal + return propose(rng, proposal, model) + end +end +function propose(rng::Random.AbstractRNG, proposals::AbstractArray{<:Proposal}, model::DensityModel, ts) + return map(proposals, ts) do proposal, t + return propose(rng, proposal, model, t) + end +end + +@generated function propose(rng::Random.AbstractRNG, proposals::NamedTuple{names}, model::DensityModel) where {names} + isempty(names) && return :(NamedTuple()) + expr = Expr(:tuple) + expr.args = Any[:($name = propose(rng, proposals.$name, model)) for name in names] + return expr +end + +@generated function propose( + rng::Random.AbstractRNG, + proposals::NamedTuple{names}, + model::DensityModel, + ts, +) where {names} + isempty(names) && return :(NamedTuple()) + expr = Expr(:tuple) + expr.args = Any[:($name = propose(rng, proposals.$name, model, ts.$name)) for name in names] + return expr +end + +# ------------------------------------------------------------------------------------------ +# logratio_proposal_density methods +# +# This is logic for the fully general case (non-symmetric, i.e. breaking detailed balance). +# For the Proposals currently implemented, these methods will be overridden by ones that +# simply return 0. + +function logratio_proposal_density(proposal::Proposal, state, candidate) + return q(proposal, state, candidate) - q(proposal, candidate, state) +end + +function logratio_proposal_density(sampler::MetropolisHastings, transition_prev::AbstractTransition, candidate) + return logratio_proposal_density(sampler.proposal, transition_prev.params, candidate) +end + +# type stable implementation for `NamedTuple`s +function logratio_proposal_density( + proposals::NamedTuple{names}, + states::NamedTuple, + candidates::NamedTuple, +) where {names} + if @generated + args = map(names) do name + :(logratio_proposal_density( + proposals[$(QuoteNode(name))], + states[$(QuoteNode(name))], + candidates[$(QuoteNode(name))], + )) + end + return :(+($(args...))) + else + return sum(names) do name + return logratio_proposal_density(proposals[name], states[name], candidates[name]) + end + end +end + +# use recursion for `Tuple`s to ensure type stability +logratio_proposal_density(proposals::Tuple{}, states::Tuple, candidates::Tuple) = 0 + +function logratio_proposal_density(proposals::Tuple{<:Proposal}, states::Tuple, candidates::Tuple) + return logratio_proposal_density(first(proposals), first(states), first(candidates)) +end + +function logratio_proposal_density(proposals::Tuple, states::Tuple, candidates::Tuple) + valfirst = logratio_proposal_density(first(proposals), first(states), first(candidates)) + valtail = logratio_proposal_density(Base.tail(proposals), Base.tail(states), Base.tail(candidates)) + return valfirst + valtail +end + +# fallback for general iterators (arrays etc.) - possibly not type stable! +function logratio_proposal_density(proposals, states, candidates) + return sum(zip(proposals, states, candidates)) do (proposal, state, candidate) + return logratio_proposal_density(proposal, state, candidate) + end +end + +# ------------------------------------------------------------------------------------------ +# bundle_samples() methods +# +# Used to capture sampler history, using MCMCChains. +# See AdvancedMH.jl/src/mcmcchains-connect.jl + +""" +A basic MCMC.Chains constructor that works with the `AbstractTransition` struct defined +in AdvancedMH.jl. +""" +function AbstractMCMC.bundle_samples( + ts::Vector{<:AbstractTransition}, + model::DensityModel, + sampler::MHSampler, + state, + chain_type::Type{Vector{NamedTuple}}; + param_names = missing, + kwargs..., +) + # Check if we received any parameter names. + if ismissing(param_names) + param_names = ["param_$i" for i in 1:length(keys(ts[1].params))] + else + # Deepcopy to be thread safe. + param_names = deepcopy(param_names) + end + + push!(param_names, "lp") + + # Turn all the transitions into a vector-of-NamedTuple. + ks = tuple(Symbol.(param_names)...) + nts = [NamedTuple{ks}(tuple(t.params..., t.lp)) for t in ts] + return nts +end + +function AbstractMCMC.bundle_samples( + ts::Vector{<:AbstractTransition}, + model::DensityModel, + sampler::MHSampler, + state, + chain_type::Type{MCMCChains.Chains}; + discard_initial = 0, + thinning = 1, + param_names = missing, + kwargs..., +) + # Turn all the transitions into a vector-of-vectors. + vals = [vcat(t.params, t.lp) for t in ts] + + # Check if we received any parameter names. + if ismissing(param_names) + param_names = [Symbol(:param_, i) for i in 1:length(keys(ts[1].params))] + else + # Generate new array to be thread safe. + param_names = Symbol.(param_names) + end + + # Bundle everything up and return a Chains struct. + return MCMCChains.Chains( + vals, + vcat(param_names, [:lp]), + (parameters = param_names, internals = [:lp]); + start = discard_initial + 1, + thin = thinning, + ) +end + +function AbstractMCMC.bundle_samples( + ts::Vector{<:Transition{<:NamedTuple}}, + model::DensityModel, + sampler::MHSampler, + state, + chain_type::Type{MCMCChains.Chains}; + discard_initial = 0, + thinning = 1, + param_names = missing, + kwargs..., +) + # Convert to a Vector{NamedTuple} first + nts = + AbstractMCMC.bundle_samples(ts, model, sampler, state, Vector{NamedTuple}; param_names = param_names, kwargs...) + + # Get all the keys + all_keys = unique(mapreduce(collect ∘ keys, vcat, nts)) + + # Push linearized draws onto array + trygetproperty(thing, key) = key in keys(thing) ? getproperty(thing, key) : missing + vals = map(nt -> [trygetproperty(nt, k) for k in all_keys], nts) + + # Check if we received any parameter names. + if ismissing(param_names) + param_names = all_keys + else + # Generate new array to be thread safe. + param_names = Symbol.(param_names) + end + + # Bundle everything up and return a Chains struct. + return MCMCChains.Chains( + vals, + param_names, + (parameters = param_names, internals = [:lp]); + start = discard_initial + 1, + thin = thinning, + ) +end + +""" +If the element type of `ts`` is `NamedTuple`s, just use the names in the struct. +""" +function AbstractMCMC.bundle_samples( + ts::Vector{<:Transition{<:NamedTuple}}, + model::DensityModel, + sampler::MHSampler, + state, + chain_type::Type{Vector{NamedTuple}}; + param_names = missing, + kwargs..., +) + # Extract NamedTuples + nts = map(x -> merge(x.params, (lp = x.lp,)), ts) + return nts +end diff --git a/src/MarkovChainMonteCarlo.jl b/src/MarkovChainMonteCarlo.jl index 36d40e141..ba1aa6bdb 100644 --- a/src/MarkovChainMonteCarlo.jl +++ b/src/MarkovChainMonteCarlo.jl @@ -22,52 +22,45 @@ abstract type AbstractMCMCAlgo end struct RandomWalkMetropolis <: AbstractMCMCAlgo end """ - MCMC{FT<:AbstractFloat, IT<:Int} +$(DocStringExtensions.TYPEDEF) -Structure to organize MCMC parameters and data +Structure to organize MCMC parameters and data. # Fields -$(DocStringExtensions.FIELDS) +$(DocStringExtensions.TYPEDFIELDS) """ -mutable struct MCMC{FT<:AbstractFloat, IT<:Int} - "a single sample from the observations. Can e.g. be picked from an Observation struct using get_obs_sample" - obs_sample::Vector{FT} - "covariance of the observational noise" +mutable struct MCMC{FT <: AbstractFloat, IT <: Int} + "A single sample from the observations. Can e.g. be picked from an `Obs` struct using `get_obs_sample`." + obs_sample::AbstractVector{FT} + "Covariance of the observational noise." obs_noise_cov::Union{AbstractMatrix{FT}, UniformScaling{FT}} - "array of length N_parameters with the parameters' prior distributions" + "Array of length *N\\_parameters* with the parameters' prior distributions." prior::ParameterDistribution - "MCMC step size" + "MCMC step size." step::FT - "Number of MCMC steps that are considered burnin" + "Number of MCMC steps that are considered burnin." burnin::IT - "the current parameters" + "The current parameters." param::AbstractVector{FT} - "Array of accepted MCMC parameter samples. The histogram of these samples gives an approximation of the posterior distribution of the parameters. param_dim x n_samples" + "Array of accepted MCMC parameter samples (*param\\_dim* × *n\\_samples*). The histogram of these samples gives an approximation of the posterior distribution of the parameters." posterior::AbstractMatrix{FT} - "the (current) value of the logarithm of the posterior (= log_likelihood + log_prior of the current parameters)" + "The current value of the logarithm of the posterior (= `log_likelihood` + `log_prior` of the current parameters)." log_posterior::Union{FT, Nothing} - "iteration/step of the MCMC" + "Iteration/step of the MCMC." iter::IT - "number of accepted proposals" + "Number of accepted proposals." accept::IT - "MCMC algorithm to use - currently implemented: 'rwm' (random walk Metropolis), 'pCN' (preconditioned Crank-Nicholson)" + "MCMC algorithm to use. Currently implemented: `'rmw'` (random walk Metropolis), `'pCN'` (preconditioned Crank-Nicholson)." algtype::String "Random number generator object (algorithm + seed) used for sampling and noise, for reproducibility." rng::Random.AbstractRNG end """ - MCMC(obs_sample::Vector{FT}, - obs_noise_cov::Array{FT, 2}, - priors::Array{Prior, 1}, - step::FT, - param_init::Vector{FT}, - max_iter::IT, - algtype::String, - burnin::IT, - -where max_iter is the number of MCMC steps to perform (e.g., 100_000) +$(DocStringExtensions.TYPEDSIGNATURES) +Constructor for [`MCMC`](@ref). +- `max_iter` - The number of MCMC steps to perform (e.g., 100_000). """ function MCMC( obs_sample::AbstractVector{FT}, @@ -82,17 +75,17 @@ function MCMC( standardize = false, norm_factor::Union{AbstractVector{FT}, Nothing} = nothing, truncate_svd = 1.0, - rng::Random.AbstractRNG = Random.GLOBAL_RNG -) where {FT<:AbstractFloat, IT<:Int} + rng::Random.AbstractRNG = Random.GLOBAL_RNG, +) where {FT <: AbstractFloat, IT <: Int} param_init_copy = deepcopy(param_init) # Standardize MCMC input? println(obs_sample) println(obs_noise_cov) if standardize - obs_sample = obs_sample ./ norm_factor; - cov_norm_factor = norm_factor .* norm_factor; - obs_noise_cov = obs_noise_cov ./ cov_norm_factor; + obs_sample = obs_sample ./ norm_factor + cov_norm_factor = norm_factor .* norm_factor + obs_noise_cov = obs_noise_cov ./ cov_norm_factor end println(obs_sample) println(obs_noise_cov) @@ -100,7 +93,7 @@ function MCMC( # We need to transform obs_sample into the correct space if svdflag println("Applying SVD to decorrelating outputs, if not required set svdflag=false") - obs_sample, _ = Emulators.svd_transform(obs_sample, obs_noise_cov; truncate_svd=truncate_svd) + obs_sample, _ = Emulators.svd_transform(obs_sample, obs_noise_cov; truncate_svd = truncate_svd) else println("Assuming independent outputs.") end @@ -114,11 +107,14 @@ function MCMC( iter = 1 accept = 0 if !(algtype in ("rwm", "pCN")) - error("Unrecognized method: ", algtype, - "Currently implemented methods: 'rwm' = random walk metropolis, ", - "'pCN' = preconditioned Crank-Nicholson") + error( + "Unrecognized method: ", + algtype, + "Currently implemented methods: 'rwm' = random walk metropolis, ", + "'pCN' = preconditioned Crank-Nicholson", + ) end - MCMC{FT,IT}( + MCMC{FT, IT}( obs_sample, obs_noise_cov, prior, @@ -130,18 +126,18 @@ function MCMC( iter, accept, algtype, - rng + rng, ) end -function reset_with_step!(mcmc::MCMC{FT,IT}, step::FT) where {FT<:AbstractFloat, IT<:Int} +function reset_with_step!(mcmc::MCMC{FT, IT}, step::FT) where {FT <: AbstractFloat, IT <: Int} # reset to beginning with new stepsize mcmc.step = step mcmc.log_posterior = nothing mcmc.iter = 1 mcmc.accept = 0 - mcmc.posterior[:,2:end] = zeros(size(mcmc.posterior[:,2:end])) + mcmc.posterior[:, 2:end] = zeros(size(mcmc.posterior[:, 2:end])) mcmc.param[:] = mcmc.posterior[:, 1] end @@ -149,7 +145,7 @@ end function get_posterior(mcmc::MCMC) #Return a parameter distributions object parameter_slices = batch(mcmc.prior) - posterior_samples = [Samples(mcmc.posterior[slice,mcmc.burnin+1:end]) for slice in parameter_slices] + posterior_samples = [Samples(mcmc.posterior[slice, (mcmc.burnin + 1):end]) for slice in parameter_slices] flattened_constraints = get_all_constraints(mcmc.prior) parameter_constraints = [flattened_constraints[slice] for slice in parameter_slices] #live in same space as prior parameter_names = get_name(mcmc.prior) #the same parameters as in prior @@ -159,10 +155,10 @@ function get_posterior(mcmc::MCMC) end function mcmc_sample!( - mcmc::MCMC{FT,IT}, - g::AbstractVector{FT}, - gcov::Union{AbstractMatrix{FT}, UniformScaling{FT}} -) where {FT<:AbstractFloat, IT<:Int} + mcmc::MCMC{FT, IT}, + g::AbstractVector{FT}, + gcov::Union{AbstractMatrix{FT}, UniformScaling{FT}}, +) where {FT <: AbstractFloat, IT <: Int} if mcmc.algtype == "rwm" log_posterior = log_likelihood(mcmc, g, gcov) + log_prior(mcmc) elseif mcmc.algtype == "pCN" @@ -181,11 +177,11 @@ function mcmc_sample!( p_accept = exp(log_posterior - mcmc.log_posterior) if p_accept > rand(mcmc.rng, Distributions.Uniform(0, 1)) - mcmc.posterior[:,1 + mcmc.iter] = mcmc.param + mcmc.posterior[:, 1 + mcmc.iter] = mcmc.param mcmc.log_posterior = log_posterior mcmc.accept = mcmc.accept + 1 else - mcmc.posterior[:, 1 + mcmc.iter[1]] = mcmc.posterior[:,mcmc.iter] + mcmc.posterior[:, 1 + mcmc.iter[1]] = mcmc.posterior[:, mcmc.iter] end mcmc.param = proposal(mcmc)[:] mcmc.iter = mcmc.iter + 1 @@ -193,37 +189,37 @@ function mcmc_sample!( end function mcmc_sample!( - mcmc::MCMC{FT,IT}, + mcmc::MCMC{FT, IT}, g::AbstractVector{FT}, - gvar::AbstractVector{FT} -) where {FT<:AbstractFloat, IT<:Int} + gvar::AbstractVector{FT}, +) where {FT <: AbstractFloat, IT <: Int} return mcmc_sample!(mcmc, g, Diagonal(gvar)) end -function accept_ratio(mcmc::MCMC{FT, IT}) where {FT<:AbstractFloat, IT<:Int} +function accept_ratio(mcmc::MCMC{FT, IT}) where {FT <: AbstractFloat, IT <: Int} return convert(FT, mcmc.accept) / mcmc.iter end function log_likelihood( - mcmc::MCMC{FT,IT}, + mcmc::MCMC{FT, IT}, g::AbstractVector{FT}, - gcov::Union{AbstractMatrix{FT}, UniformScaling{FT}} -) where {FT<:AbstractFloat, IT<:Int} + gcov::Union{AbstractMatrix{FT}, UniformScaling{FT}}, +) where {FT <: AbstractFloat, IT <: Int} log_rho = 0.0 #if gcov == nothing # diff = g - mcmc.obs_sample # log_rho[1] = -FT(0.5) * diff' * (mcmc.obs_noise_cov \ diff) #else - # det(log(Γ)) - # Ill-posed numerically for ill-conditioned covariance matrices with det≈0 - #log_gpfidelity = -FT(0.5) * log(det(Diagonal(gvar))) # = -0.5 * sum(log.(gvar)) - # Well-posed numerically for ill-conditioned covariance matrices with det≈0 - #full_cov = Diagonal(gvar) - eigs = eigvals(gcov) - log_gpfidelity = -FT(0.5) * sum(log.(eigs)) + # det(log(Γ)) + # Ill-posed numerically for ill-conditioned covariance matrices with det≈0 + #log_gpfidelity = -FT(0.5) * log(det(Diagonal(gvar))) # = -0.5 * sum(log.(gvar)) + # Well-posed numerically for ill-conditioned covariance matrices with det≈0 + #full_cov = Diagonal(gvar) + eigs = eigvals(gcov) + log_gpfidelity = -FT(0.5) * sum(log.(eigs)) # Combine got log_rho - diff = g - mcmc.obs_sample + diff = g - mcmc.obs_sample log_rho = -FT(0.5) * diff' * (gcov \ diff) + log_gpfidelity #end return log_rho @@ -231,7 +227,7 @@ end function log_prior(mcmc::MCMC) - return get_logpdf(mcmc.prior,mcmc.param) + return get_logpdf(mcmc.prior, mcmc.param) end @@ -240,12 +236,12 @@ function proposal(mcmc::MCMC) prop_dist = MvNormal(zeros(length(mcmc.param)), proposal_covariance) if mcmc.algtype == "rwm" - sample = mcmc.posterior[:,1 + mcmc.iter] .+ mcmc.step * rand(mcmc.rng, prop_dist) + sample = mcmc.posterior[:, 1 + mcmc.iter] .+ mcmc.step * rand(mcmc.rng, prop_dist) elseif mcmc.algtype == "pCN" # Use prescription in Beskos et al (2017) "Geometric MCMC for infinite-dimensional # inverse problems." for relating ρ to Euler stepsize: - ρ = (1 - mcmc.step/4) / (1 + mcmc.step/4) - sample = ρ * mcmc.posterior[:,1 + mcmc.iter] .+ sqrt(1 - ρ^2) * rand(mcmc.rng, prop_dist) + ρ = (1 - mcmc.step / 4) / (1 + mcmc.step / 4) + sample = ρ * mcmc.posterior[:, 1 + mcmc.iter] .+ sqrt(1 - ρ^2) * rand(mcmc.rng, prop_dist) else error("Unrecognized algtype: ", mcmc.algtype) end @@ -253,7 +249,11 @@ function proposal(mcmc::MCMC) end -function find_mcmc_step!(mcmc_test::MCMC{FT,IT}, em::Emulator{FT}; max_iter::IT=2000) where {FT<:AbstractFloat, IT<:Int} +function find_mcmc_step!( + mcmc_test::MCMC{FT, IT}, + em::Emulator{FT}; + max_iter::IT = 2000, +) where {FT <: AbstractFloat, IT <: Int} step = mcmc_test.step mcmc_accept = false doubled = false @@ -268,7 +268,7 @@ function find_mcmc_step!(mcmc_test::MCMC{FT,IT}, em::Emulator{FT}; max_iter::IT= while mcmc_accept == false param = reshape(mcmc_test.param, :, 1) - em_pred, em_predvar = predict(em, param ) + em_pred, em_predvar = predict(em, param) if ndims(em_predvar[1]) != 0 mcmc_sample!(mcmc_test, vec(em_pred), diag(em_predvar[1])) else @@ -278,12 +278,10 @@ function find_mcmc_step!(mcmc_test::MCMC{FT,IT}, em::Emulator{FT}; max_iter::IT= if it % max_iter == 0 countmcmc += 1 acc_ratio = accept_ratio(mcmc_test) - println("iteration ", it, "; acceptance rate = ", acc_ratio, - ", current parameters ", param) + println("iteration ", it, "; acceptance rate = ", acc_ratio, ", current parameters ", param) flush(stdout) if countmcmc == 20 - println("failed to choose suitable stepsize in ", countmcmc, - "iterations") + println("failed to choose suitable stepsize in ", countmcmc, "iterations") exit() end it = 0 @@ -296,7 +294,7 @@ function find_mcmc_step!(mcmc_test::MCMC{FT,IT}, em::Emulator{FT}; max_iter::IT= step *= 0.5 reset_with_step!(mcmc_test, step) halved = true - elseif acc_ratio>0.35 + elseif acc_ratio > 0.35 step *= 2.0 reset_with_step!(mcmc_test, step) doubled = true @@ -315,11 +313,7 @@ function find_mcmc_step!(mcmc_test::MCMC{FT,IT}, em::Emulator{FT}; max_iter::IT= end -function sample_posterior!( - mcmc::MCMC{FT, IT}, - em::Emulator{FT}, - max_iter::IT -) where {FT<:AbstractFloat, IT<:Int} +function sample_posterior!(mcmc::MCMC{FT, IT}, em::Emulator{FT}, max_iter::IT) where {FT <: AbstractFloat, IT <: Int} for mcmcit in 1:max_iter param = reshape(mcmc.param, :, 1) # test predictions (param is 1 x N_parameters) diff --git a/src/Utilities.jl b/src/Utilities.jl index 2f98f1777..eb354842e 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -1,5 +1,6 @@ module Utilities +using DocStringExtensions using LinearAlgebra using Statistics using StatsBase @@ -15,22 +16,22 @@ export orig2zscore export zscore2orig """ - extract_training_points(ekp::EnsembleKalmanProcess{FT, IT, P}, N_train_iterations::IT) where {FT,IT, P} +$(DocStringExtensions.TYPEDSIGNATURES) -Extract the training points needed to train the Gaussian Process Regression. +Extract the training points needed to train the Gaussian process regression. - - `ekp` - EnsembleKalmanProcess holding the parameters and the data that were produced - during the Ensemble Kalman (EK) process - - `N_train_iterations` - Number of EK layers/iterations to train on +- `ekp` - EnsembleKalmanProcess holding the parameters and the data that were produced + during the Ensemble Kalman (EK) process. +- `N_train_iterations` - Number of EK layers/iterations to train on. """ function get_training_points(ekp::EnsembleKalmanProcess{FT, IT, P}, N_train_iterations::IT) where {FT, IT, P} # Note u[end] does not have an equivalent g - iter_range = get_N_iterations(ekp) - N_train_iterations + 1 : get_N_iterations(ekp) + iter_range = (get_N_iterations(ekp) - N_train_iterations + 1):get_N_iterations(ekp) - u_tp=[] - g_tp=[] + u_tp = [] + g_tp = [] for i in iter_range push!(u_tp, get_u(ekp, i)) #N_parameters x N_ens push!(g_tp, get_g(ekp, i)) #N_data x N_ens @@ -38,16 +39,16 @@ function get_training_points(ekp::EnsembleKalmanProcess{FT, IT, P}, N_train_iter u_tp = hcat(u_tp...) # N_parameters x (N_ek_it x N_ensemble)] g_tp = hcat(g_tp...) # N_data x (N_ek_it x N_ensemble) - training_points = PairedDataContainer(u_tp,g_tp,data_are_columns=true) - + training_points = PairedDataContainer(u_tp, g_tp, data_are_columns = true) + return training_points end """ - get_obs_sample(rng::AbstractRNG, obs::Observation; rng_seed::Union{IT, Nothing} = nothing) +$(DocStringExtensions.TYPEDSIGNATURES) -Return one random sample from the observations (for use in the MCMC) +Return a random sample from the observations, for use in the MCMC. - `rng` - optional RNG object used to pick random sample; defaults to `Random.GLOBAL_RNG`. - `obs` - Observation struct with the observations (extract will pick one @@ -55,8 +56,8 @@ Return one random sample from the observations (for use in the MCMC) - `rng_seed` - optional kwarg; if provided, used to re-seed `rng` before sampling. """ function get_obs_sample( - rng::Random.AbstractRNG, - obs::Observation; + rng::Random.AbstractRNG, + obs::Observation; rng_seed::Union{IT, Nothing} = nothing, ) where {IT <: Int} # Ensuring reproducibility of the sampled parameter values: @@ -64,15 +65,13 @@ function get_obs_sample( if rng_seed !== nothing rng = Random.seed!(rng, rng_seed) end - row_idxs = StatsBase.sample(rng, axes(obs.samples, 1), 1; replace=false, ordered=false) + row_idxs = StatsBase.sample(rng, axes(obs.samples, 1), 1; replace = false, ordered = false) return obs.samples[row_idxs...] end # first arg optional; defaults to GLOBAL_RNG (as in Random, StatsBase) get_obs_sample(obs::Observation; kwargs...) = get_obs_sample(Random.GLOBAL_RNG, obs; kwargs...) -function orig2zscore(X::AbstractVector{FT}, - mean::AbstractVector{FT}, - std::AbstractVector{FT}) where {FT} +function orig2zscore(X::AbstractVector{FT}, mean::AbstractVector{FT}, std::AbstractVector{FT}) where {FT} # Compute the z scores of a vector X using the given mean # and std Z = zeros(size(X)) @@ -82,40 +81,34 @@ function orig2zscore(X::AbstractVector{FT}, return Z end -function orig2zscore(X::AbstractMatrix{FT}, - mean::AbstractVector{FT}, - std::AbstractVector{FT}) where {FT} +function orig2zscore(X::AbstractMatrix{FT}, mean::AbstractVector{FT}, std::AbstractVector{FT}) where {FT} # Compute the z scores of matrix X using the given mean and # std. Transformation is applied column-wise. Z = zeros(size(X)) n_cols = size(X)[2] for i in 1:n_cols - Z[:,i] = (X[:,i] .- mean[i]) ./ std[i] + Z[:, i] = (X[:, i] .- mean[i]) ./ std[i] end return Z end -function zscore2orig(Z::AbstractVector{FT}, - mean::AbstractVector{FT}, - std::AbstractVector{FT}) where {FT} +function zscore2orig(Z::AbstractVector{FT}, mean::AbstractVector{FT}, std::AbstractVector{FT}) where {FT} # Transform X (a vector of z scores) back to the original # values X = zeros(size(Z)) for i in 1:length(X) - X[i] = Z[i] .* std[i] .+ mean[i] + X[i] = Z[i] .* std[i] .+ mean[i] end return X end -function zscore2orig(Z::AbstractMatrix{FT}, - mean::AbstractVector{FT}, - std::AbstractVector{FT}) where {FT} +function zscore2orig(Z::AbstractMatrix{FT}, mean::AbstractVector{FT}, std::AbstractVector{FT}) where {FT} X = zeros(size(Z)) # Transform X (a matrix of z scores) back to the original # values. Transformation is applied column-wise. n_cols = size(Z)[2] for i in 1:n_cols - X[:,i] = Z[:,i] .* std[i] .+ mean[i] + X[:, i] = Z[:, i] .* std[i] .+ mean[i] end return X end diff --git a/test/Emulator/runtests.jl b/test/Emulator/runtests.jl index ef35e5dee..8e57f383c 100644 --- a/test/Emulator/runtests.jl +++ b/test/Emulator/runtests.jl @@ -1,7 +1,7 @@ # Import modules using Random using Test -using Statistics +using Statistics using Distributions using LinearAlgebra @@ -14,11 +14,12 @@ struct MLTester <: Emulators.MachineLearningTool end function constructor_tests( iopairs::PairedDataContainer, Σ::Union{AbstractMatrix{FT}, UniformScaling{FT}, Nothing}, - norm_factors, decomposition -) where {FT<:AbstractFloat} + norm_factors, + decomposition, +) where {FT <: AbstractFloat} # [2.] test Normalization - input_mean = vec(mean(get_inputs(iopairs), dims=2)) #column vector - sqrt_inv_input_cov = sqrt(inv(Symmetric(cov(get_inputs(iopairs), dims=2)))) + input_mean = vec(mean(get_inputs(iopairs), dims = 2)) #column vector + sqrt_inv_input_cov = sqrt(inv(Symmetric(cov(get_inputs(iopairs), dims = 2)))) norm_inputs = Emulators.normalize(get_inputs(iopairs), input_mean, sqrt_inv_input_cov) @test norm_inputs == sqrt_inv_input_cov * (get_inputs(iopairs) .- input_mean) @@ -27,21 +28,23 @@ function constructor_tests( @test_throws ErrorException emulator = Emulator( mlt, iopairs, - obs_noise_cov=Σ, - normalize_inputs=true, - standardize_outputs=false, - truncate_svd=1.0) + obs_noise_cov = Σ, + normalize_inputs = true, + standardize_outputs = false, + truncate_svd = 1.0, + ) #build a known type, with defaults gp = GaussianProcess(GPJL()) emulator = Emulator( gp, iopairs, - obs_noise_cov=Σ, - normalize_inputs=false, - standardize_outputs=false, - truncate_svd=1.0) - + obs_noise_cov = Σ, + normalize_inputs = false, + standardize_outputs = false, + truncate_svd = 1.0, + ) + # compare SVD/norm/stand with stored emulator version test_decomp = emulator.decomposition @test test_decomp.V == decomposition.V #(use [:,:] to make it an array) @@ -52,10 +55,11 @@ function constructor_tests( emulator2 = Emulator( gp, iopairs, - obs_noise_cov=Σ, - normalize_inputs=true, - standardize_outputs=false, - truncate_svd=1.0) + obs_noise_cov = Σ, + normalize_inputs = true, + standardize_outputs = false, + truncate_svd = 1.0, + ) train_inputs = get_inputs(emulator2.training_pairs) @test train_inputs == norm_inputs train_inputs2 = Emulators.normalize(emulator2, get_inputs(iopairs)) @@ -66,74 +70,80 @@ function constructor_tests( emulator3 = Emulator( gp, iopairs, - obs_noise_cov=Σ, - normalize_inputs=false, - standardize_outputs=true, - standardize_outputs_factors=norm_factors, - truncate_svd=1.0) + obs_noise_cov = Σ, + normalize_inputs = false, + standardize_outputs = true, + standardize_outputs_factors = norm_factors, + truncate_svd = 1.0, + ) #standardized and decorrelated (sd) data s_y, s_y_cov = Emulators.standardize(get_outputs(iopairs), Σ, norm_factors) sd_train_outputs = get_outputs(emulator3.training_pairs) - sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(emulator3.decomposition.S)) + sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(emulator3.decomposition.S)) decorrelated_s_y = sqrt_singular_values_inv * emulator3.decomposition.Vt * s_y @test decorrelated_s_y == sd_train_outputs end @testset "Emulators" begin #build some quick data + noise - m=50 - d=6 + m = 50 + d = 6 x = rand(3, m) #R^3 y = rand(d, m) #R^5 - + # "noise" - μ = zeros(d) - Σ = rand(d,d) - Σ = Σ'*Σ + μ = zeros(d) + Σ = rand(d, d) + Σ = Σ' * Σ noise_samples = rand(MvNormal(μ, Σ), m) y += noise_samples - - iopairs = PairedDataContainer(x,y,data_are_columns=true) + + iopairs = PairedDataContainer(x, y, data_are_columns = true) @test get_inputs(iopairs) == x @test get_outputs(iopairs) == y - + # [1.] test SVD (2D y version) test_SVD = svd(Σ) - transformed_y, decomposition = Emulators.svd_transform(y, Σ, truncate_svd=1.0) - @test_throws MethodError Emulators.svd_transform(y, Σ[:,1], truncate_svd=1.0) - @test test_SVD.V[:,:] == decomposition.V #(use [:,:] to make it an array) + transformed_y, decomposition = Emulators.svd_transform(y, Σ, truncate_svd = 1.0) + @test_throws MethodError Emulators.svd_transform(y, Σ[:, 1], truncate_svd = 1.0) + @test test_SVD.V[:, :] == decomposition.V #(use [:,:] to make it an array) @test test_SVD.Vt == decomposition.Vt @test test_SVD.S == decomposition.S @test size(transformed_y) == size(y) # 1D y version - transformed_y, decomposition2 = Emulators.svd_transform(y[:, 1], Σ, truncate_svd=1.0) - @test_throws MethodError Emulators.svd_transform(y[:,1], Σ[:,1], truncate_svd=1.0) + transformed_y, decomposition2 = Emulators.svd_transform(y[:, 1], Σ, truncate_svd = 1.0) + @test_throws MethodError Emulators.svd_transform(y[:, 1], Σ[:, 1], truncate_svd = 1.0) @test size(transformed_y) == size(y[:, 1]) - + # Reverse SVD - y_new, y_cov_new = Emulators.svd_reverse_transform_mean_cov(reshape(transformed_y, (d,1)), ones(d,1), decomposition2) - @test_throws MethodError Emulators.svd_reverse_transform_mean_cov(transformed_y, ones(d,1), decomposition2) - @test_throws MethodError Emulators.svd_reverse_transform_mean_cov(reshape(transformed_y,(d,1)), ones(d), decomposition2) + y_new, y_cov_new = + Emulators.svd_reverse_transform_mean_cov(reshape(transformed_y, (d, 1)), ones(d, 1), decomposition2) + @test_throws MethodError Emulators.svd_reverse_transform_mean_cov(transformed_y, ones(d, 1), decomposition2) + @test_throws MethodError Emulators.svd_reverse_transform_mean_cov( + reshape(transformed_y, (d, 1)), + ones(d), + decomposition2, + ) @test size(y_new)[1] == size(y[:, 1])[1] - @test y_new ≈ y[:,1] + @test y_new ≈ y[:, 1] @test y_cov_new[1] ≈ Σ - + # Truncation - transformed_y, trunc_decomposition = Emulators.svd_transform(y[:, 1], Σ, truncate_svd=0.95) + transformed_y, trunc_decomposition = Emulators.svd_transform(y[:, 1], Σ, truncate_svd = 0.95) trunc_size = size(trunc_decomposition.S)[1] @test test_SVD.S[1:trunc_size] == trunc_decomposition.S @test size(transformed_y)[1] == trunc_size - + # [3.] test Standardization norm_factors = 10.0 - norm_factors = fill(norm_factors, size(y[:,1])) # must be size of output dim + norm_factors = fill(norm_factors, size(y[:, 1])) # must be size of output dim + + s_y, s_y_cov = Emulators.standardize(get_outputs(iopairs), Σ, norm_factors) + @test s_y == get_outputs(iopairs) ./ norm_factors + @test s_y_cov == Σ ./ (norm_factors .* norm_factors') - s_y, s_y_cov = Emulators.standardize(get_outputs(iopairs),Σ,norm_factors) - @test s_y == get_outputs(iopairs)./norm_factors - @test s_y_cov == Σ ./ (norm_factors.*norm_factors') - @testset "Emulators - dense Array Σ" begin constructor_tests(iopairs, Σ, norm_factors, decomposition) @@ -142,14 +152,15 @@ end emulator4 = Emulator( gp, iopairs, - obs_noise_cov=Σ, - normalize_inputs=false, - standardize_outputs=false, - truncate_svd=0.9) + obs_noise_cov = Σ, + normalize_inputs = false, + standardize_outputs = false, + truncate_svd = 0.9, + ) trunc_size = size(emulator4.decomposition.S)[1] @test test_SVD.S[1:trunc_size] == emulator4.decomposition.S end - + @testset "Emulators - Diagonal Σ" begin x = rand(3, m) #R^3 y = rand(d, m) #R^5 @@ -157,11 +168,11 @@ end Σ = Diagonal(rand(d) .+ 1.0) noise_samples = rand(MvNormal(zeros(d), Σ), m) y += noise_samples - iopairs = PairedDataContainer(x, y, data_are_columns=true) - _, decomposition = Emulators.svd_transform(y, Σ, truncate_svd=1.0) + iopairs = PairedDataContainer(x, y, data_are_columns = true) + _, decomposition = Emulators.svd_transform(y, Σ, truncate_svd = 1.0) constructor_tests(iopairs, Σ, norm_factors, decomposition) - end + end @testset "Emulators - UniformScaling Σ" begin x = rand(3, m) #R^3 @@ -170,9 +181,9 @@ end Σ = UniformScaling(1.3) noise_samples = rand(MvNormal(zeros(d), Σ), m) y += noise_samples - iopairs = PairedDataContainer(x, y, data_are_columns=true) - _, decomposition = Emulators.svd_transform(y, Σ, truncate_svd=1.0) + iopairs = PairedDataContainer(x, y, data_are_columns = true) + _, decomposition = Emulators.svd_transform(y, Σ, truncate_svd = 1.0) constructor_tests(iopairs, Σ, norm_factors, decomposition) - end + end end diff --git a/test/GaussianProcess/runtests.jl b/test/GaussianProcess/runtests.jl index f7acf0dbc..3907a28e1 100644 --- a/test/GaussianProcess/runtests.jl +++ b/test/GaussianProcess/runtests.jl @@ -2,12 +2,12 @@ using Random using Test using GaussianProcesses -using Statistics +using Statistics using Distributions using ScikitLearn using LinearAlgebra -@sk_import gaussian_process : GaussianProcessRegressor -@sk_import gaussian_process.kernels : (RBF, WhiteKernel, ConstantKernel) +@sk_import gaussian_process:GaussianProcessRegressor +@sk_import (gaussian_process.kernels):(RBF, WhiteKernel, ConstantKernel) using CalibrateEmulateSample.Emulators using CalibrateEmulateSample.DataContainers @@ -27,15 +27,15 @@ using CalibrateEmulateSample.DataContainers x = reshape(2.0 * π * rand(n), 1, n) # predictors/features: 1 x n y = reshape(sin.(x) + 0.05 * randn(n)', 1, n) # predictands/targets: 1 x n - iopairs = PairedDataContainer(x,y,data_are_columns=true) + iopairs = PairedDataContainer(x, y, data_are_columns = true) # Construct kernel: # Squared exponential kernel (note that hyperparameters are on log scale) # with observational noise GPkernel = SE(log(1.0), log(1.0)) # These will be the test inputs at which predictions are made - new_inputs = reshape([0.0, π/2, π, 3*π/2, 2*π], 1, 5) - + new_inputs = reshape([0.0, π / 2, π, 3 * π / 2, 2 * π], 1, 5) + # Fit Gaussian Process Regression models. # GaussianProcesses.jl (GPJL) provides two predict functions, predict_y # (which predicts the random variable y(θ)) and predict_y (which predicts @@ -46,92 +46,83 @@ using CalibrateEmulateSample.DataContainers ## GaussianProcess 1: GPJL, predict_y unnormalized gppackage = GPJL() pred_type = YType() - - gp1 = GaussianProcess( - gppackage; - kernel=GPkernel, - noise_learn=true, - prediction_type=pred_type) + + gp1 = GaussianProcess(gppackage; kernel = GPkernel, noise_learn = true, prediction_type = pred_type) @test gp1.kernel == GPkernel @test gp1.noise_learn == true @test gp1.prediction_type == pred_type - + em1 = Emulator( gp1, iopairs, - obs_noise_cov=nothing, - normalize_inputs=false, - standardize_outputs=false, - truncate_svd=1.0) + obs_noise_cov = nothing, + normalize_inputs = false, + standardize_outputs = false, + truncate_svd = 1.0, + ) Emulators.optimize_hyperparameters!(em1) μ1, σ1² = Emulators.predict(em1, new_inputs) - - @test vec(μ1) ≈ [0.0, 1.0, 0.0, -1.0, 0.0] atol=0.3 - @test size(μ1)==(1,5) - @test vec(σ1²) ≈ [0.017, 0.003, 0.004, 0.004, 0.009] atol=1e-2 + + @test vec(μ1) ≈ [0.0, 1.0, 0.0, -1.0, 0.0] atol = 0.3 + @test size(μ1) == (1, 5) + @test vec(σ1²) ≈ [0.017, 0.003, 0.004, 0.004, 0.009] atol = 1e-2 # GaussianProcess 2: GPJL, predict_f pred_type = FType() - gp2 = GaussianProcess( - gppackage; - kernel=GPkernel, - noise_learn=true, - prediction_type=pred_type) + gp2 = GaussianProcess(gppackage; kernel = GPkernel, noise_learn = true, prediction_type = pred_type) em2 = Emulator( gp2, iopairs, - obs_noise_cov=nothing, - normalize_inputs=false, - standardize_outputs=false, - truncate_svd=1.0) + obs_noise_cov = nothing, + normalize_inputs = false, + standardize_outputs = false, + truncate_svd = 1.0, + ) Emulators.optimize_hyperparameters!(em2) μ2, σ2² = Emulators.predict(em2, new_inputs) # predict_y and predict_f should give the same mean - @test μ2 ≈ μ1 atol=1e-6 + @test μ2 ≈ μ1 atol = 1e-6 # GaussianProcess 3: SKLJL gppackage = SKLJL() pred_type = YType() - var = ConstantKernel(constant_value=1.0) + var = ConstantKernel(constant_value = 1.0) se = RBF(1.0) GPkernel = var * se - gp3 = GaussianProcess( - gppackage; - kernel=GPkernel, - noise_learn=true, - prediction_type=pred_type) + gp3 = GaussianProcess(gppackage; kernel = GPkernel, noise_learn = true, prediction_type = pred_type) em3 = Emulator( gp3, iopairs, - obs_noise_cov=nothing, - normalize_inputs=false, - standardize_outputs=false, - truncate_svd=1.0) + obs_noise_cov = nothing, + normalize_inputs = false, + standardize_outputs = false, + truncate_svd = 1.0, + ) Emulators.optimize_hyperparameters!(em3) - + #gp3 = GaussianProcess(iopairs, gppackage; GPkernel=GPkernel, obs_noise_cov=nothing, # normalized=false, noise_learn=true, # truncate_svd=1.0, standardize=false, # prediction_type=pred_type, norm_factor=nothing) - + μ3, σ3² = Emulators.predict(em3, new_inputs) - @test vec(μ3) ≈ [0.0, 1.0, 0.0, -1.0, 0.0] atol=0.3 - @test vec(σ3²) ≈ [0.016, 0.002, 0.003, 0.004, 0.003] atol=1e-2 - + @test vec(μ3) ≈ [0.0, 1.0, 0.0, -1.0, 0.0] atol = 0.3 + @test vec(σ3²) ≈ [0.016, 0.002, 0.003, 0.004, 0.003] atol = 1e-2 + # ------------------------------------------------------------------------- # Test case 2: 2D input, 2D output @@ -142,58 +133,55 @@ using CalibrateEmulateSample.DataContainers # Generate training data m = 80 # number of training points - + p = 2 # input dim d = 2 # output dim X = 2.0 * π * rand(p, m) - + # G(x1, x2) g1x = sin.(X[1, :]) .+ cos.(X[2, :]) g2x = sin.(X[1, :]) .- cos.(X[2, :]) - gx = zeros(2,m) - gx[1,:] = g1x - gx[2,:] = g2x + gx = zeros(2, m) + gx[1, :] = g1x + gx[2, :] = g2x # Add noise η - μ = zeros(d) + μ = zeros(d) Σ = 0.1 * [[0.8, 0.2] [0.2, 0.5]] # d x d noise_samples = rand(MvNormal(μ, Σ), m) - + # y = G(x) + η Y = gx .+ noise_samples - iopairs2 = PairedDataContainer(X,Y,data_are_columns=true) + iopairs2 = PairedDataContainer(X, Y, data_are_columns = true) @test get_inputs(iopairs2) == X @test get_outputs(iopairs2) == Y - gp4 = GaussianProcess( - gppackage; - kernel=nothing, - noise_learn=true, - prediction_type=pred_type) + gp4 = GaussianProcess(gppackage; kernel = nothing, noise_learn = true, prediction_type = pred_type) em4 = Emulator( gp4, iopairs2, - obs_noise_cov=Σ, - normalize_inputs=true, - standardize_outputs=false, - truncate_svd=1.0) + obs_noise_cov = Σ, + normalize_inputs = true, + standardize_outputs = false, + truncate_svd = 1.0, + ) Emulators.optimize_hyperparameters!(em4) - + new_inputs = zeros(2, 4) - new_inputs[:, 2] = [π/2, π] - new_inputs[:, 3] = [π, π/2] - new_inputs[:, 4] = [3*π/2, 2*π] - - μ4, σ4² = Emulators.predict(em4, new_inputs, transform_to_real=true) - - @test μ4[:, 1] ≈ [1.0, -1.0] atol=0.25 - @test μ4[:, 2] ≈ [0.0, 2.0] atol=0.25 - @test μ4[:, 3] ≈ [0.0, 0.0] atol=0.25 - @test μ4[:, 4] ≈ [0.0, -2.0] atol=0.25 - @test length(σ4²) == size(new_inputs,2) - @test size(σ4²[1]) == (d, d) + new_inputs[:, 2] = [π / 2, π] + new_inputs[:, 3] = [π, π / 2] + new_inputs[:, 4] = [3 * π / 2, 2 * π] + + μ4, σ4² = Emulators.predict(em4, new_inputs, transform_to_real = true) + + @test μ4[:, 1] ≈ [1.0, -1.0] atol = 0.25 + @test μ4[:, 2] ≈ [0.0, 2.0] atol = 0.25 + @test μ4[:, 3] ≈ [0.0, 0.0] atol = 0.25 + @test μ4[:, 4] ≈ [0.0, -2.0] atol = 0.25 + @test length(σ4²) == size(new_inputs, 2) + @test size(σ4²[1]) == (d, d) end diff --git a/test/MarkovChainMonteCarlo/runtests.jl b/test/MarkovChainMonteCarlo/runtests.jl index 141ba599d..74bce5818 100644 --- a/test/MarkovChainMonteCarlo/runtests.jl +++ b/test/MarkovChainMonteCarlo/runtests.jl @@ -16,10 +16,10 @@ function test_gp_data(; rng_seed = 41, n = 20, var_y = 0.05, rest...) # that's tested in test/GaussianProcesses/runtests.jl Case 1 n = 40 # number of training points x = 2π * rand(rng, Float64, (1, n)) # predictors/features: 1 × n - σ2_y = reshape([var_y],1,1) + σ2_y = reshape([var_y], 1, 1) y = sin.(x) + rand(rng, Normal(0, σ2_y[1]), (1, n)) # predictands/targets: 1 × n - return y, σ2_y, PairedDataContainer(x, y, data_are_columns=true), rng + return y, σ2_y, PairedDataContainer(x, y, data_are_columns = true), rng end function test_gp_prior() @@ -27,7 +27,7 @@ function test_gp_prior() umin = -1.0 umax = 6.0 #prior = [Priors.Prior(Uniform(umin, umax), "u")] # prior on u - prior_dist = Parameterized(Normal(0,1)) + prior_dist = Parameterized(Normal(0, 1)) prior_constraint = bounded(umin, umax) prior_name = "u" return ParameterDistribution(prior_dist, prior_constraint, prior_name) @@ -40,14 +40,15 @@ function test_gp_1(y, σ2_y, iopairs::PairedDataContainer; norm_factor = nothing # Squared exponential kernel (note that hyperparameters are on log scale) # with observational noise GPkernel = SE(log(1.0), log(1.0)) - gp = GaussianProcess( - gppackage; - kernel=GPkernel, noise_learn=true, prediction_type=pred_type - ) + gp = GaussianProcess(gppackage; kernel = GPkernel, noise_learn = true, prediction_type = pred_type) em = Emulator( - gp, iopairs; - obs_noise_cov=σ2_y, normalize_inputs=false, - standardize_outputs=false, standardize_outputs_factors=norm_factor, truncate_svd=1.0 + gp, + iopairs; + obs_noise_cov = σ2_y, + normalize_inputs = false, + standardize_outputs = false, + standardize_outputs_factors = norm_factor, + truncate_svd = 1.0, ) Emulators.optimize_hyperparameters!(em) return em @@ -60,22 +61,30 @@ function test_gp_2(y, σ2_y, iopairs::PairedDataContainer; norm_factor = nothing # Squared exponential kernel (note that hyperparameters are on log scale) # with observational noise GPkernel = SE(log(1.0), log(1.0)) - gp = GaussianProcess( - gppackage; - kernel=GPkernel, noise_learn=true, prediction_type=pred_type - ) + gp = GaussianProcess(gppackage; kernel = GPkernel, noise_learn = true, prediction_type = pred_type) em = Emulator( - gp, iopairs; - obs_noise_cov=σ2_y, normalize_inputs=false, - standardize_outputs=true, standardize_outputs_factors=norm_factor, truncate_svd=0.9, + gp, + iopairs; + obs_noise_cov = σ2_y, + normalize_inputs = false, + standardize_outputs = true, + standardize_outputs_factors = norm_factor, + truncate_svd = 0.9, ) Emulators.optimize_hyperparameters!(em) return em end -function mcmc_gp_test_step(prior:: ParameterDistribution, σ2_y, em::Emulator; - mcmc_alg = "rwm", obs_sample = 1.0, param_init = 3.0, step = 0.5, norm_factor=nothing, - rng=Random.GLOBAL_RNG +function mcmc_gp_test_step( + prior::ParameterDistribution, + σ2_y, + em::Emulator; + mcmc_alg = "rwm", + obs_sample = 1.0, + param_init = 3.0, + step = 0.5, + norm_factor = nothing, + rng = Random.GLOBAL_RNG, ) obs_sample = reshape(collect(obs_sample), 1) # scalar or Vector -> Vector param_init = reshape(collect(param_init), 1) # scalar or Vector -> Vector @@ -83,18 +92,35 @@ function mcmc_gp_test_step(prior:: ParameterDistribution, σ2_y, em::Emulator; burnin = 0 step = 0.5 # first guess max_iter = 5000 - mcmc_test = MCMC(obs_sample, σ2_y, prior, step, param_init, max_iter, mcmc_alg, burnin; - svdflag=true, standardize=false, truncate_svd=1.0, - norm_factor=nothing, rng=rng - ) + mcmc_test = MCMC( + obs_sample, + σ2_y, + prior, + step, + param_init, + max_iter, + mcmc_alg, + burnin; + svdflag = true, + standardize = false, + truncate_svd = 1.0, + norm_factor = nothing, + rng = rng, + ) new_step = find_mcmc_step!(mcmc_test, em) return new_step end function mcmc_gp_test_template( - prior:: ParameterDistribution, σ2_y, em::Emulator; - mcmc_alg = "rwm", obs_sample = 1.0, param_init = 3.0, step = 0.5, norm_factor=nothing, - rng=Random.GLOBAL_RNG + prior::ParameterDistribution, + σ2_y, + em::Emulator; + mcmc_alg = "rwm", + obs_sample = 1.0, + param_init = 3.0, + step = 0.5, + norm_factor = nothing, + rng = Random.GLOBAL_RNG, ) obs_sample = reshape(collect(obs_sample), 1) # scalar or Vector -> Vector param_init = reshape(collect(param_init), 1) # scalar or Vector -> Vector @@ -102,32 +128,48 @@ function mcmc_gp_test_template( burnin = 1000 max_iter = 100_000 # Now begin the actual MCMC - mcmc = MCMC(obs_sample, σ2_y, prior, step, param_init, max_iter, mcmc_alg, burnin; - svdflag=true, standardize=false, truncate_svd=1.0, norm_factor=norm_factor, - rng=rng - ) + mcmc = MCMC( + obs_sample, + σ2_y, + prior, + step, + param_init, + max_iter, + mcmc_alg, + burnin; + svdflag = true, + standardize = false, + truncate_svd = 1.0, + norm_factor = norm_factor, + rng = rng, + ) sample_posterior!(mcmc, em, max_iter) - posterior_distribution = get_posterior(mcmc) + posterior_distribution = get_posterior(mcmc) #post_mean = mean(posterior, dims=1)[1] posterior_mean = mean(posterior_distribution) # We had svdflag=true, so the MCMC stores a transformed sample, # 1.0/sqrt(0.05) * obs_sample ≈ 4.472 - @test mcmc.obs_sample ≈ (obs_sample ./ sqrt(σ2_y[1,1])) atol=1e-2 + @test mcmc.obs_sample ≈ (obs_sample ./ sqrt(σ2_y[1, 1])) atol = 1e-2 @test mcmc.obs_noise_cov == σ2_y @test mcmc.burnin == burnin @test mcmc.algtype == mcmc_alg @test mcmc.iter[1] == max_iter + 1 @test get_n_samples(posterior_distribution)[prior.names[1]] == max_iter - burnin + 1 @test ndims(posterior_distribution) == length(param_init) - @test_throws Exception MCMC(obs_sample, σ2_y, prior, step, param_init, - max_iter, "gibbs", burnin) + @test_throws Exception MCMC(obs_sample, σ2_y, prior, step, param_init, max_iter, "gibbs", burnin) return posterior_mean[1] end function mcmc_gp_standardization_test_template( - prior:: ParameterDistribution, σ2_y, em::Emulator; - mcmc_alg = "rwm", obs_sample = 1.0, param_init = 3.0, step = 0.5, norm_factor=nothing, - rng=Random.GLOBAL_RNG + prior::ParameterDistribution, + σ2_y, + em::Emulator; + mcmc_alg = "rwm", + obs_sample = 1.0, + param_init = 3.0, + step = 0.5, + norm_factor = nothing, + rng = Random.GLOBAL_RNG, ) obs_sample = reshape(collect(obs_sample), 1) # scalar or Vector -> Vector param_init = reshape(collect(param_init), 1) # scalar or Vector -> Vector @@ -135,15 +177,37 @@ function mcmc_gp_standardization_test_template( burnin = 1000 max_iter = 100_000 # Standardization and truncation - mcmc_test = MCMC(obs_sample, σ2_y, prior, step, param_init, max_iter, mcmc_alg, burnin; - svdflag=true, standardize=true, norm_factor=norm_factor, - truncate_svd=0.9, rng=rng - ) + mcmc_test = MCMC( + obs_sample, + σ2_y, + prior, + step, + param_init, + max_iter, + mcmc_alg, + burnin; + svdflag = true, + standardize = true, + norm_factor = norm_factor, + truncate_svd = 0.9, + rng = rng, + ) # Now begin the actual MCMC - mcmc = MCMC(obs_sample, σ2_y, prior, step, param_init, max_iter, mcmc_alg, burnin; - svdflag=true, standardize=false, truncate_svd=1.0, norm_factor=norm_factor, - rng=rng - ) + mcmc = MCMC( + obs_sample, + σ2_y, + prior, + step, + param_init, + max_iter, + mcmc_alg, + burnin; + svdflag = true, + standardize = false, + truncate_svd = 1.0, + norm_factor = norm_factor, + rng = rng, + ) sample_posterior!(mcmc, em, max_iter) posterior_mean = mean(get_posterior(mcmc)) return posterior_mean[1] @@ -154,47 +218,59 @@ end prior = test_gp_prior() y, σ2_y, iopairs, rng = test_gp_data(; rng_seed = 456, n = 20, var_y = 0.05) em = test_gp_1(y, σ2_y, iopairs) - mcmc_params = Dict(:mcmc_alg => "rwm", :obs_sample => 1.0, :param_init => 3.0, - :step => 0.5, :norm_factor=> nothing, :rng => rng) + mcmc_params = Dict( + :mcmc_alg => "rwm", + :obs_sample => 1.0, + :param_init => 3.0, + :step => 0.5, + :norm_factor => nothing, + :rng => rng, + ) new_step = mcmc_gp_test_step(prior, σ2_y, em; mcmc_params...) - @test isapprox(new_step, 1.0; atol=0.1) + @test isapprox(new_step, 1.0; atol = 0.1) posterior_mean_1 = mcmc_gp_test_template(prior, σ2_y, em; mcmc_params...) # difference between mean_1 and ground truth comes from MCMC convergence and GP sampling - @test isapprox(posterior_mean_1, π/2; atol=4e-1) + @test isapprox(posterior_mean_1, π / 2; atol = 4e-1) norm_factor = 10.0 - norm_factor = fill(norm_factor, size(y[:,1])) # must be size of output dim - em = test_gp_2(y, σ2_y, iopairs; norm_factor=norm_factor) - posterior_mean_2 = mcmc_gp_standardization_test_template(prior, σ2_y, em; - mcmc_params..., norm_factor=norm_factor) + norm_factor = fill(norm_factor, size(y[:, 1])) # must be size of output dim + em = test_gp_2(y, σ2_y, iopairs; norm_factor = norm_factor) + posterior_mean_2 = + mcmc_gp_standardization_test_template(prior, σ2_y, em; mcmc_params..., norm_factor = norm_factor) # difference between mean_1 and mean_2 only from MCMC convergence - @test isapprox(posterior_mean_2, posterior_mean_1; atol=0.1) + @test isapprox(posterior_mean_2, posterior_mean_1; atol = 0.1) end @testset "GP & pCN" begin prior = test_gp_prior() y, σ2_y, iopairs, rng = test_gp_data(; rng_seed = 456, n = 20, var_y = 0.05) em = test_gp_1(y, σ2_y, iopairs) - mcmc_params = Dict(:mcmc_alg => "pCN", :obs_sample => 1.0, :param_init => 3.0, - :step => 0.5, :norm_factor=> nothing, :rng => rng) + mcmc_params = Dict( + :mcmc_alg => "pCN", + :obs_sample => 1.0, + :param_init => 3.0, + :step => 0.5, + :norm_factor => nothing, + :rng => rng, + ) new_step = mcmc_gp_test_step(prior, σ2_y, em; mcmc_params...) - @test isapprox(new_step, 1.0; atol=0.1) + @test isapprox(new_step, 1.0; atol = 0.1) posterior_mean_1 = mcmc_gp_test_template(prior, σ2_y, em; mcmc_params...) # difference between mean_1 and ground truth comes from MCMC convergence and GP sampling - @test isapprox(posterior_mean_1, π/2; atol=4e-1) + @test isapprox(posterior_mean_1, π / 2; atol = 4e-1) norm_factor = 10.0 - norm_factor = fill(norm_factor, size(y[:,1])) # must be size of output dim - em = test_gp_2(y, σ2_y, iopairs; norm_factor=norm_factor) - posterior_mean_2 = mcmc_gp_standardization_test_template(prior, σ2_y, em; - mcmc_params..., norm_factor=norm_factor) + norm_factor = fill(norm_factor, size(y[:, 1])) # must be size of output dim + em = test_gp_2(y, σ2_y, iopairs; norm_factor = norm_factor) + posterior_mean_2 = + mcmc_gp_standardization_test_template(prior, σ2_y, em; mcmc_params..., norm_factor = norm_factor) # difference between mean_1 and mean_2 only from MCMC convergence - @test isapprox(posterior_mean_2, posterior_mean_1; atol=0.1) + @test isapprox(posterior_mean_2, posterior_mean_1; atol = 0.1) end end diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index 1d6ec93fd..bec9283ef 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -13,46 +13,48 @@ using CalibrateEmulateSample.DataContainers # Seed for pseudo-random number generator rng = Random.MersenneTwister(41) - arr = vcat([i*ones(3)' for i in 1:5]...) - arr_t = permutedims(arr,(2,1)) + arr = vcat([i * ones(3)' for i in 1:5]...) + arr_t = permutedims(arr, (2, 1)) data_names = ["d1", "d2", "d3"] obs = Observation(arr_t, data_names) #data must be columns as default sample = get_obs_sample(rng, obs) @test sample == [5.0, 5.0, 5.0] - mean_arr = dropdims(mean(arr, dims=1), dims=1) - std_arr = dropdims(std(arr, dims=1), dims=1) + mean_arr = dropdims(mean(arr, dims = 1), dims = 1) + std_arr = dropdims(std(arr, dims = 1), dims = 1) @test mean_arr == [3.0, 3.0, 3.0] - @test std_arr ≈ [1.58, 1.58, 1.58] atol=1e-2 + @test std_arr ≈ [1.58, 1.58, 1.58] atol = 1e-2 z_arr = orig2zscore(arr, mean_arr, std_arr) - z_arr_test = [-1.265 -1.265 -1.265; - -0.632 -0.632 -0.632; - 0.0 0.0 0.0; - 0.632 0.632 0.632; - 1.265 1.265 1.265] - @test z_arr ≈ z_arr_test atol=1e-2 + z_arr_test = [ + -1.265 -1.265 -1.265 + -0.632 -0.632 -0.632 + 0.0 0.0 0.0 + 0.632 0.632 0.632 + 1.265 1.265 1.265 + ] + @test z_arr ≈ z_arr_test atol = 1e-2 orig_arr = zscore2orig(z_arr, mean_arr, std_arr) - @test orig_arr ≈ arr atol=1e-5 + @test orig_arr ≈ arr atol = 1e-5 v = vec([1.0 2.0 3.0]) mean_v = vec([0.5 1.5 2.5]) std_v = vec([0.5 0.5 0.5]) z_v = orig2zscore(v, mean_v, std_v) println(z_v) - @test z_v ≈ [1.0, 1.0, 1.0] atol=1e-5 + @test z_v ≈ [1.0, 1.0, 1.0] atol = 1e-5 orig_v = zscore2orig(z_v, mean_v, std_v) - @test orig_v ≈ v atol=1e-5 + @test orig_v ≈ v atol = 1e-5 # test get_training_points # first create the EnsembleKalmanProcess n_ens = 10 dim_obs = 3 dim_par = 2 - initial_ensemble = randn(rng, dim_par,n_ens)#params are cols + initial_ensemble = randn(rng, dim_par, n_ens)#params are cols y_obs = randn(rng, dim_obs) - Γy = Matrix{Float64}(I,dim_obs,dim_obs) - ekp = EnsembleKalmanProcesses.EnsembleKalmanProcess(initial_ensemble, y_obs, Γy, Inversion(), rng=rng) - g_ens = randn(rng, dim_obs,n_ens) # data are cols + Γy = Matrix{Float64}(I, dim_obs, dim_obs) + ekp = EnsembleKalmanProcesses.EnsembleKalmanProcess(initial_ensemble, y_obs, Γy, Inversion(), rng = rng) + g_ens = randn(rng, dim_obs, n_ens) # data are cols EnsembleKalmanProcesses.update_ensemble!(ekp, g_ens) training_points = get_training_points(ekp, 1) @test get_inputs(training_points) ≈ initial_ensemble diff --git a/test/runtests.jl b/test/runtests.jl index 17dae84ab..c48b2aa9d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -19,18 +19,14 @@ end function has_submodule(sm) any(ARGS) do a a == sm && return true - first(split(a, '/')) == sm && return true - return false + first(split(a, '/')) == sm && return true + return false end end - for submodule in ["Emulator", - "GaussianProcess", - "MarkovChainMonteCarlo", - "Utilities"] - + for submodule in ["Emulator", "GaussianProcess", "MarkovChainMonteCarlo", "Utilities"] if all_tests || has_submodule(submodule) || "CalibrateEmulateSample" in ARGS include_test(submodule) end - end + end end