diff --git a/examples/Cloudy/Cloudy_example.jl b/examples/Cloudy/Cloudy_example.jl
index 94adb7d2e..798772d37 100644
--- a/examples/Cloudy/Cloudy_example.jl
+++ b/examples/Cloudy/Cloudy_example.jl
@@ -206,7 +206,8 @@ kern2 = Mat52Ard(len2, 0.0)
 white = Noise(log(2.0))
 GPkernel =  kern1 + kern2 + white
 # Get training points
-u_tp, g_tp = Utilities.extract_GP_tp(ekiobj, N_iter)
+#u_tp, g_tp = Utilities.extract_GP_tp(ekiobj, N_iter)
+u_tp, g_tp = Utilities.get_training_points(ekiobj, N_iter)
 normalized = true
 gpobj = GaussianProcessEmulator.GaussianProcess(u_tp, g_tp, gppackage; GPkernel=GPkernel, 
                          obs_noise_cov=Γy, normalized=normalized, 
diff --git a/examples/Cloudy/Manifest.toml b/examples/Cloudy/Manifest.toml
deleted file mode 100644
index 0b63d996d..000000000
--- a/examples/Cloudy/Manifest.toml
+++ /dev/null
@@ -1,1264 +0,0 @@
-# This file is machine-generated - editing it directly is not advised
-
-[[AbstractFFTs]]
-deps = ["LinearAlgebra"]
-git-tree-sha1 = "051c95d6836228d120f5f4b984dd5aba1624f716"
-uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
-version = "0.5.0"
-
-[[AbstractTrees]]
-deps = ["Markdown"]
-git-tree-sha1 = "33e450545eaf7699da1a6e755f9ea65f14077a45"
-uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
-version = "0.3.3"
-
-[[Adapt]]
-deps = ["LinearAlgebra"]
-git-tree-sha1 = "0fac443759fa829ed8066db6cf1077d888bb6573"
-uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "2.0.2"
-
-[[ArnoldiMethod]]
-deps = ["DelimitedFiles", "LinearAlgebra", "Random", "SparseArrays", "StaticArrays", "Test"]
-git-tree-sha1 = "2b6845cea546604fb4dca4e31414a6a59d39ddcd"
-uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
-version = "0.0.4"
-
-[[Arpack]]
-deps = ["Arpack_jll", "Libdl", "LinearAlgebra"]
-git-tree-sha1 = "2ff92b71ba1747c5fdd541f8fc87736d82f40ec9"
-uuid = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
-version = "0.4.0"
-
-[[Arpack_jll]]
-deps = ["Libdl", "OpenBLAS_jll", "Pkg"]
-git-tree-sha1 = "e214a9b9bd1b4e1b4f15b22c0994862b66af7ff7"
-uuid = "68821587-b530-5797-8361-c406ea357684"
-version = "3.5.0+3"
-
-[[ArrayInterface]]
-deps = ["LinearAlgebra", "Requires", "SparseArrays"]
-git-tree-sha1 = "066d1e7a9eb4873660791db7f0d8c7902600b81c"
-uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "2.11.0"
-
-[[ArrayLayouts]]
-deps = ["FillArrays", "LinearAlgebra"]
-git-tree-sha1 = "6f6e33efac70fc24c1f2a654a090b7af01690ffe"
-uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
-version = "0.3.7"
-
-[[AxisAlgorithms]]
-deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"]
-git-tree-sha1 = "a4d07a1c313392a77042855df46c5f534076fab9"
-uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950"
-version = "1.0.0"
-
-[[BandedMatrices]]
-deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "Random", "SparseArrays"]
-git-tree-sha1 = "eaf98fa821ab26c5825fa3a054db735755c335da"
-uuid = "aae01518-5342-5314-be14-df237901396f"
-version = "0.15.15"
-
-[[Base64]]
-uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
-
-[[BinDeps]]
-deps = ["Libdl", "Pkg", "SHA", "URIParser", "Unicode"]
-git-tree-sha1 = "46cf2c1668ad07aba5a9d331bdeea994a1f13856"
-uuid = "9e28174c-4ba2-5203-b857-d8d62c4213ee"
-version = "1.0.1"
-
-[[BinaryProvider]]
-deps = ["Libdl", "Logging", "SHA"]
-git-tree-sha1 = "ecdec412a9abc8db54c0efc5548c64dfce072058"
-uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
-version = "0.5.10"
-
-[[BoundaryValueDiffEq]]
-deps = ["BandedMatrices", "DiffEqBase", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "NLsolve", "Reexport", "SparseArrays"]
-git-tree-sha1 = "3f2969de608af70db755cee9d4490a7294a6afc3"
-uuid = "764a87c0-6b3e-53db-9096-fe964310641d"
-version = "2.5.0"
-
-[[Bzip2_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "3663bfffede2ef41358b6fc2e1d8a6d50b3c3904"
-uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
-version = "1.0.6+2"
-
-[[CEnum]]
-git-tree-sha1 = "215a9aa4a1f23fbd05b92769fdd62559488d70e9"
-uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
-version = "0.4.1"
-
-[[CalibrateEmulateSample]]
-deps = ["Cloudy", "Conda", "DifferentialEquations", "Distributions", "DocStringExtensions", "GaussianProcesses", "LinearAlgebra", "Optim", "Pkg", "Plots", "PyCall", "Random", "ScikitLearn", "Statistics", "StatsBase", "StatsPlots", "Sundials", "Test"]
-git-tree-sha1 = "8afef38a7a09401c9bc5d7f3cf952c05afe071e5"
-repo-rev = "master"
-repo-url = "https://github.com/CliMA/CalibrateEmulateSample.jl.git"
-uuid = "95e48a1f-0bec-4818-9538-3db4340308e3"
-version = "0.1.0"
-
-[[CanonicalTraits]]
-deps = ["MLStyle"]
-git-tree-sha1 = "1728d2633a206e931daf818309d90eb39c1ce93b"
-uuid = "a603d957-0e48-4f86-8fbd-0b7bc66df689"
-version = "0.2.1"
-
-[[CategoricalArrays]]
-deps = ["DataAPI", "Future", "JSON", "Missings", "Printf", "Statistics", "Unicode"]
-git-tree-sha1 = "a6c17353ee38ddab30e73dcfaa1107752de724ec"
-uuid = "324d7699-5711-5eae-9e2f-1d82baa6b597"
-version = "0.8.1"
-
-[[ChainRulesCore]]
-deps = ["MuladdMacro"]
-git-tree-sha1 = "971b03f25bdf2acab79f1c51afc717f9dccf43c2"
-uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-version = "0.9.5"
-
-[[Cloudy]]
-deps = ["Coverage", "DifferentialEquations", "DocStringExtensions", "ForwardDiff", "HCubature", "LinearAlgebra", "Optim", "PyPlot", "SpecialFunctions", "TaylorSeries", "Test"]
-git-tree-sha1 = "7cb3ce76f5a63d17eb358b6030f7664931c75857"
-uuid = "9e3b23bb-e7cc-4b94-886c-65de2234ba87"
-version = "0.1.0"
-
-[[Clustering]]
-deps = ["Distances", "LinearAlgebra", "NearestNeighbors", "Printf", "SparseArrays", "Statistics", "StatsBase"]
-git-tree-sha1 = "b11c8d607af357776a046889a7c32567d05f1319"
-uuid = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5"
-version = "0.14.1"
-
-[[ColorSchemes]]
-deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random", "StaticArrays"]
-git-tree-sha1 = "7a15e3690529fd1042f0ab954dff7445b1efc8a5"
-uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
-version = "3.9.0"
-
-[[ColorTypes]]
-deps = ["FixedPointNumbers", "Random"]
-git-tree-sha1 = "6e7aa35d0294f647bb9c985ccc34d4f5d371a533"
-uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
-version = "0.10.6"
-
-[[Colors]]
-deps = ["ColorTypes", "FixedPointNumbers", "InteractiveUtils", "Reexport"]
-git-tree-sha1 = "5639e44833cfcf78c6a73fbceb4da75611d312cd"
-uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
-version = "0.12.3"
-
-[[Combinatorics]]
-git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860"
-uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
-version = "1.0.2"
-
-[[CommonSubexpressions]]
-deps = ["MacroTools", "Test"]
-git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
-uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
-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 = "a6a8197ae253f2c1a22b2ae17c2dfaf5812c03aa"
-uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "3.13.0"
-
-[[CompilerSupportLibraries_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "7c4f882c41faa72118841185afc58a2eb00ef612"
-uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "0.3.3+0"
-
-[[Conda]]
-deps = ["JSON", "VersionParsing"]
-git-tree-sha1 = "7a58bb32ce5d85f8bf7559aa7c2842f9aecf52fc"
-uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d"
-version = "1.4.1"
-
-[[ConsoleProgressMonitor]]
-deps = ["Logging", "ProgressMeter"]
-git-tree-sha1 = "3ab7b2136722890b9af903859afcf457fa3059e8"
-uuid = "88cd18e8-d9cc-4ea6-8889-5259c0d15c8b"
-version = "0.1.2"
-
-[[ConstructionBase]]
-git-tree-sha1 = "a2a6a5fea4d6f730ec4c18a76d27ec10e8ec1c50"
-uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
-version = "1.0.0"
-
-[[Contour]]
-deps = ["StaticArrays"]
-git-tree-sha1 = "81685fee51fc5168898e3cbd8b0f01506cd9148e"
-uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
-version = "0.5.4"
-
-[[Coverage]]
-deps = ["CoverageTools", "HTTP", "JSON", "LibGit2", "MbedTLS"]
-git-tree-sha1 = "aad69c0c2d8080427f507c36b289c77fe465c9fe"
-uuid = "a2441757-f6aa-5fb2-8edb-039e3f45d037"
-version = "1.1.1"
-
-[[CoverageTools]]
-deps = ["LibGit2"]
-git-tree-sha1 = "328154f6062b0053588295df6bffe22e1436597a"
-uuid = "c36e975a-824b-4404-a568-ef97ca766997"
-version = "1.1.0"
-
-[[CpuId]]
-deps = ["Markdown", "Test"]
-git-tree-sha1 = "f0464e499ab9973b43c20f8216d088b61fda80c6"
-uuid = "adafc99b-e345-5852-983c-f28acb93d879"
-version = "0.2.2"
-
-[[DataAPI]]
-git-tree-sha1 = "176e23402d80e7743fc26c19c681bfb11246af32"
-uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
-version = "1.3.0"
-
-[[DataFrames]]
-deps = ["CategoricalArrays", "Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "Missings", "PooledArrays", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
-git-tree-sha1 = "d4436b646615928b634b37e99a3288588072f851"
-uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-version = "0.21.4"
-
-[[DataStructures]]
-deps = ["InteractiveUtils", "OrderedCollections"]
-git-tree-sha1 = "edad9434967fdc0a2631a65d902228400642120c"
-uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-version = "0.17.19"
-
-[[DataValueInterfaces]]
-git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
-uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
-version = "1.0.0"
-
-[[DataValues]]
-deps = ["DataValueInterfaces", "Dates"]
-git-tree-sha1 = "d88a19299eba280a6d062e135a43f00323ae70bf"
-uuid = "e7dc6d0d-1eca-5fa6-8ad6-5aecde8b7ea5"
-version = "0.4.13"
-
-[[Dates]]
-deps = ["Printf"]
-uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
-
-[[DelayDiffEq]]
-deps = ["DataStructures", "DiffEqBase", "LinearAlgebra", "Logging", "OrdinaryDiffEq", "Printf", "RecursiveArrayTools", "Reexport", "Roots", "UnPack"]
-git-tree-sha1 = "ede535e084589ca2ed626f4179337ac7fd610e7f"
-uuid = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb"
-version = "5.24.1"
-
-[[DelimitedFiles]]
-deps = ["Mmap"]
-uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
-
-[[DiffEqBase]]
-deps = ["ArrayInterface", "ChainRulesCore", "ConsoleProgressMonitor", "DataStructures", "Distributed", "DocStringExtensions", "FunctionWrappers", "IterativeSolvers", "IteratorInterfaceExtensions", "LabelledArrays", "LinearAlgebra", "Logging", "LoggingExtras", "MuladdMacro", "Parameters", "Printf", "ProgressLogging", "RecipesBase", "RecursiveArrayTools", "RecursiveFactorization", "Requires", "Roots", "SparseArrays", "StaticArrays", "Statistics", "SuiteSparse", "TableTraits", "TerminalLoggers", "TreeViews", "ZygoteRules"]
-git-tree-sha1 = "7f241c049129a0e765ef449b634990409286ed75"
-uuid = "2b5f629d-d688-5b77-993f-72d75c75574e"
-version = "6.40.9"
-
-[[DiffEqCallbacks]]
-deps = ["DataStructures", "DiffEqBase", "ForwardDiff", "LinearAlgebra", "NLsolve", "OrdinaryDiffEq", "RecipesBase", "RecursiveArrayTools", "StaticArrays"]
-git-tree-sha1 = "63afe1bc146559c7329b6429713916c7f9e61d2f"
-uuid = "459566f4-90b8-5000-8ac3-15dfb0a30def"
-version = "2.13.5"
-
-[[DiffEqFinancial]]
-deps = ["DiffEqBase", "DiffEqNoiseProcess", "LinearAlgebra", "Markdown", "RandomNumbers"]
-git-tree-sha1 = "db08e0def560f204167c58fd0637298e13f58f73"
-uuid = "5a0ffddc-d203-54b0-88ba-2c03c0fc2e67"
-version = "2.4.0"
-
-[[DiffEqJump]]
-deps = ["ArrayInterface", "Compat", "DataStructures", "DiffEqBase", "FunctionWrappers", "LinearAlgebra", "Parameters", "PoissonRandom", "Random", "RandomNumbers", "RecursiveArrayTools", "StaticArrays", "Statistics", "TreeViews"]
-git-tree-sha1 = "c9baaba9b1ee1407473a2daac1ca0ffdcb72253d"
-uuid = "c894b116-72e5-5b58-be3c-e6d8d4ac2b12"
-version = "6.9.3"
-
-[[DiffEqNoiseProcess]]
-deps = ["DataStructures", "DiffEqBase", "Distributions", "LinearAlgebra", "PoissonRandom", "Random", "RandomNumbers", "RecipesBase", "RecursiveArrayTools", "Requires", "ResettableStacks", "StaticArrays", "Statistics"]
-git-tree-sha1 = "52ede62cc53338d1b46c2278ef85837dc90ead62"
-uuid = "77a26b50-5914-5dd7-bc55-306e6241c503"
-version = "5.1.1"
-
-[[DiffEqPhysics]]
-deps = ["DiffEqBase", "DiffEqCallbacks", "ForwardDiff", "LinearAlgebra", "Printf", "Random", "RecipesBase", "RecursiveArrayTools", "Reexport", "StaticArrays"]
-git-tree-sha1 = "9f4b98590c63bf202d12ae20783203baad3dd27a"
-uuid = "055956cb-9e8b-5191-98cc-73ae4a59e68a"
-version = "3.2.0"
-
-[[DiffResults]]
-deps = ["StaticArrays"]
-git-tree-sha1 = "da24935df8e0c6cf28de340b958f6aac88eaa0cc"
-uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
-version = "1.0.2"
-
-[[DiffRules]]
-deps = ["NaNMath", "Random", "SpecialFunctions"]
-git-tree-sha1 = "eb0c34204c8410888844ada5359ac8b96292cfd1"
-uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
-version = "1.0.1"
-
-[[DifferentialEquations]]
-deps = ["BoundaryValueDiffEq", "DelayDiffEq", "DiffEqBase", "DiffEqCallbacks", "DiffEqFinancial", "DiffEqJump", "DiffEqNoiseProcess", "DiffEqPhysics", "DimensionalPlotRecipes", "LinearAlgebra", "MultiScaleArrays", "OrdinaryDiffEq", "ParameterizedFunctions", "Random", "RecursiveArrayTools", "Reexport", "SteadyStateDiffEq", "StochasticDiffEq", "Sundials"]
-git-tree-sha1 = "84502233aa79c6b2d0b9d5d50d205b88a02ffd8b"
-uuid = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
-version = "6.15.0"
-
-[[DimensionalPlotRecipes]]
-deps = ["LinearAlgebra", "RecipesBase"]
-git-tree-sha1 = "af883a26bbe6e3f5f778cb4e1b81578b534c32a6"
-uuid = "c619ae07-58cd-5f6d-b883-8f17bd6a98f9"
-version = "1.2.0"
-
-[[Distances]]
-deps = ["LinearAlgebra", "Statistics"]
-git-tree-sha1 = "bed62cc5afcff16de797a9f38fb358b74071f785"
-uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
-version = "0.9.0"
-
-[[Distributed]]
-deps = ["Random", "Serialization", "Sockets"]
-uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
-
-[[Distributions]]
-deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"]
-git-tree-sha1 = "21fe8727469b4eeb5a006e0399588180f5dc9db4"
-uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
-version = "0.23.5"
-
-[[DocStringExtensions]]
-deps = ["LibGit2", "Markdown", "Pkg", "Test"]
-git-tree-sha1 = "c5714d9bcdba66389612dc4c47ed827c64112997"
-uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-version = "0.8.2"
-
-[[Documenter]]
-deps = ["Base64", "Dates", "DocStringExtensions", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
-git-tree-sha1 = "395fa1554c69735802bba37d9e7d9586fd44326c"
-uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-version = "0.24.11"
-
-[[ElasticArrays]]
-deps = ["Adapt"]
-git-tree-sha1 = "b740f7f3fa87a2fa22c93582f2dbb9110c4ece6e"
-uuid = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4"
-version = "1.2.3"
-
-[[ElasticPDMats]]
-deps = ["LinearAlgebra", "MacroTools", "PDMats", "Test"]
-git-tree-sha1 = "c6ac4406d9b4a549f7316fc746f7ccd5f1bcd2cd"
-uuid = "2904ab23-551e-5aed-883f-487f97af5226"
-version = "0.2.1"
-
-[[ExponentialUtilities]]
-deps = ["LinearAlgebra", "Printf", "Requires", "SparseArrays"]
-git-tree-sha1 = "91f7498b66205431fe3e35833cda97a22b1ab6a5"
-uuid = "d4d017d3-3776-5f7e-afef-a10c40355c18"
-version = "1.7.0"
-
-[[FFMPEG]]
-deps = ["FFMPEG_jll"]
-git-tree-sha1 = "c82bef6fc01e30d500f588cd01d29bdd44f1924e"
-uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
-version = "0.3.0"
-
-[[FFMPEG_jll]]
-deps = ["Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "LAME_jll", "LibVPX_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "Pkg", "Zlib_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
-git-tree-sha1 = "0fa07f43e5609ea54848b82b4bb330b250e9645b"
-uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
-version = "4.1.0+3"
-
-[[FFTW]]
-deps = ["AbstractFFTs", "FFTW_jll", "IntelOpenMP_jll", "Libdl", "LinearAlgebra", "MKL_jll", "Reexport"]
-git-tree-sha1 = "14536c95939aadcee44014728a459d2fe3ca9acf"
-uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
-version = "1.2.2"
-
-[[FFTW_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "6c975cd606128d45d1df432fb812d6eb10fee00b"
-uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
-version = "3.3.9+5"
-
-[[FastGaussQuadrature]]
-deps = ["LinearAlgebra", "SpecialFunctions"]
-git-tree-sha1 = "c139e3f4c75dc489a493627c7ee44befc177420f"
-uuid = "442a2c76-b920-505d-bb47-c5924d526838"
-version = "0.4.2"
-
-[[FillArrays]]
-deps = ["LinearAlgebra", "Random", "SparseArrays"]
-git-tree-sha1 = "be4180bdb27a11188d694ee3773122f4921f1a62"
-uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "0.8.13"
-
-[[FiniteDiff]]
-deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"]
-git-tree-sha1 = "43b397bf66d07d113cc2cc012dd2029d336c3894"
-uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
-version = "2.5.1"
-
-[[FixedPointNumbers]]
-deps = ["Statistics"]
-git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
-uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
-version = "0.8.4"
-
-[[Formatting]]
-deps = ["Printf"]
-git-tree-sha1 = "a0c901c29c0e7c763342751c0a94211d56c0de5c"
-uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
-version = "0.4.1"
-
-[[ForwardDiff]]
-deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "NaNMath", "Random", "SpecialFunctions", "StaticArrays"]
-git-tree-sha1 = "1d090099fb82223abc48f7ce176d3f7696ede36d"
-uuid = "f6369f11-7733-5829-9624-2563aa707210"
-version = "0.10.12"
-
-[[FreeType2_jll]]
-deps = ["Bzip2_jll", "Libdl", "Pkg", "Zlib_jll"]
-git-tree-sha1 = "7d900f32a3788d4eacac2bfa3bf5c770179c8afd"
-uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
-version = "2.10.1+2"
-
-[[FriBidi_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "2f56bee16bd0151de7b6a1eeea2ced190a2ad8d4"
-uuid = "559328eb-81f9-559d-9380-de523a88c83c"
-version = "1.0.5+3"
-
-[[FunctionWrappers]]
-git-tree-sha1 = "e4813d187be8c7b993cb7f85cbf2b7bfbaadc694"
-uuid = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
-version = "1.1.1"
-
-[[Future]]
-deps = ["Random"]
-uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
-
-[[GR]]
-deps = ["Base64", "DelimitedFiles", "HTTP", "JSON", "LinearAlgebra", "Printf", "Random", "Serialization", "Sockets", "Test", "UUIDs"]
-git-tree-sha1 = "e26c513329675092535de20cc4bb9c579c8f85a0"
-uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
-version = "0.51.0"
-
-[[GaussianProcesses]]
-deps = ["Distances", "Distributions", "Documenter", "ElasticArrays", "ElasticPDMats", "FastGaussQuadrature", "ForwardDiff", "LinearAlgebra", "Optim", "PDMats", "Printf", "ProgressMeter", "Random", "RecipesBase", "ScikitLearnBase", "SpecialFunctions", "StaticArrays", "Statistics", "StatsFuns", "Zygote"]
-git-tree-sha1 = "9b33cd8be0fd2bc59c46911def57331a7d5f6e84"
-uuid = "891a1506-143c-57d2-908e-e1f8e92e6de9"
-version = "0.12.1"
-
-[[GeneralizedGenerated]]
-deps = ["CanonicalTraits", "DataStructures", "JuliaVariables", "MLStyle"]
-git-tree-sha1 = "9a917449bebc5a241d23e13d36ea62c77129ce6e"
-uuid = "6b9d7cbe-bcb9-11e9-073f-15a7a543e2eb"
-version = "0.2.4"
-
-[[GenericSVD]]
-deps = ["LinearAlgebra"]
-git-tree-sha1 = "62909c3eda8a25b5673a367d1ad2392ebb265211"
-uuid = "01680d73-4ee2-5a08-a1aa-533608c188bb"
-version = "0.3.0"
-
-[[GeometryBasics]]
-deps = ["IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
-git-tree-sha1 = "119f32f9c2b497b49cd3f7f513b358b82660294c"
-uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
-version = "0.2.15"
-
-[[GeometryTypes]]
-deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "StaticArrays"]
-git-tree-sha1 = "34bfa994967e893ab2f17b864eec221b3521ba4d"
-uuid = "4d00f742-c7ba-57c2-abde-4428a4b178cb"
-version = "0.8.3"
-
-[[HCubature]]
-deps = ["Combinatorics", "DataStructures", "LinearAlgebra", "QuadGK", "StaticArrays"]
-git-tree-sha1 = "fc4adde0a029bf69c0da3f7c25ca24b5e0e389ed"
-uuid = "19dc6840-f33b-545b-b366-655c7e3ffd49"
-version = "1.4.0"
-
-[[HTTP]]
-deps = ["Base64", "Dates", "IniFile", "MbedTLS", "Sockets"]
-git-tree-sha1 = "eca61b35cdd8cd2fcc5eec1eda766424a995b02f"
-uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
-version = "0.8.16"
-
-[[IRTools]]
-deps = ["InteractiveUtils", "MacroTools", "Test"]
-git-tree-sha1 = "90ee39f9beaaa186e4968417ea2b8ed5673c91c0"
-uuid = "7869d1d1-7146-5819-86e3-90919afe41df"
-version = "0.3.3"
-
-[[Inflate]]
-git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
-uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
-version = "0.1.2"
-
-[[IniFile]]
-deps = ["Test"]
-git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
-uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
-version = "0.5.0"
-
-[[IntelOpenMP_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "fb8e1c7a5594ba56f9011310790e03b5384998d6"
-uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
-version = "2018.0.3+0"
-
-[[InteractiveUtils]]
-deps = ["Markdown"]
-uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
-
-[[Interpolations]]
-deps = ["AxisAlgorithms", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"]
-git-tree-sha1 = "2b7d4e9be8b74f03115e64cf36ed2f48ae83d946"
-uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
-version = "0.12.10"
-
-[[InvertedIndices]]
-deps = ["Test"]
-git-tree-sha1 = "15732c475062348b0165684ffe28e85ea8396afc"
-uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
-version = "1.0.0"
-
-[[IterTools]]
-git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
-uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
-version = "1.3.0"
-
-[[IterativeSolvers]]
-deps = ["LinearAlgebra", "Printf", "Random", "RecipesBase", "SparseArrays"]
-git-tree-sha1 = "3b7e2aac8c94444947facea7cc7ca91c49169be0"
-uuid = "42fd0dbc-a981-5370-80f2-aaf504508153"
-version = "0.8.4"
-
-[[IteratorInterfaceExtensions]]
-git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
-uuid = "82899510-4779-5014-852e-03e436cf321d"
-version = "1.0.0"
-
-[[JSON]]
-deps = ["Dates", "Mmap", "Parsers", "Unicode"]
-git-tree-sha1 = "b34d7cef7b337321e97d22242c3c2b91f476748e"
-uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
-version = "0.21.0"
-
-[[JuliaVariables]]
-deps = ["MLStyle", "NameResolution"]
-git-tree-sha1 = "8868479ff35ab98588ed0a529a9c2a4f8bda3bd6"
-uuid = "b14d175d-62b4-44ba-8fb7-3064adc8c3ec"
-version = "0.2.0"
-
-[[KernelDensity]]
-deps = ["Distributions", "FFTW", "Interpolations", "Optim", "StatsBase", "Test"]
-git-tree-sha1 = "c1048817fe5711f699abc8fabd47b1ac6ba4db04"
-uuid = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
-version = "0.5.1"
-
-[[LAME_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "221cc8998b9060677448cbb6375f00032554c4fd"
-uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d"
-version = "3.100.0+1"
-
-[[LaTeXStrings]]
-git-tree-sha1 = "de44b395389b84fd681394d4e8d39ef14e3a2ea8"
-uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
-version = "1.1.0"
-
-[[LabelledArrays]]
-deps = ["ArrayInterface", "LinearAlgebra", "MacroTools", "StaticArrays"]
-git-tree-sha1 = "5e04374019448f8509349948ab504f117e3b575a"
-uuid = "2ee39098-c373-598a-b85f-a56591580800"
-version = "1.3.0"
-
-[[Latexify]]
-deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
-git-tree-sha1 = "864527aa4d14c893fb8c51d48ef314410c88c7b9"
-uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
-version = "0.13.5"
-
-[[LeftChildRightSiblingTrees]]
-deps = ["AbstractTrees"]
-git-tree-sha1 = "71be1eb5ad19cb4f61fa8c73395c0338fd092ae0"
-uuid = "1d6d02ad-be62-4b6b-8a6d-2f90e265016e"
-version = "0.1.2"
-
-[[LibGit2]]
-uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
-
-[[LibVPX_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "e3549ca9bf35feb9d9d954f4c6a9032e92f46e7c"
-uuid = "dd192d2f-8180-539f-9fb4-cc70b1dcf69a"
-version = "1.8.1+1"
-
-[[Libdl]]
-uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
-
-[[LightGraphs]]
-deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
-git-tree-sha1 = "6f85a35d2377cb2db1bc448ed0d6340d2bb1ea64"
-uuid = "093fc24a-ae57-5d10-9952-331d41423f4d"
-version = "1.3.3"
-
-[[LineSearches]]
-deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf", "Test"]
-git-tree-sha1 = "54eb90e8dbe745d617c78dee1d6ae95c7f6f5779"
-uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
-version = "7.0.1"
-
-[[LinearAlgebra]]
-deps = ["Libdl"]
-uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-
-[[Logging]]
-uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
-
-[[LoggingExtras]]
-deps = ["Dates"]
-git-tree-sha1 = "03289aba73c0abc25ff0229bed60f2a4129cd15c"
-uuid = "e6f89c97-d47a-5376-807f-9c37f3926c36"
-version = "0.4.2"
-
-[[LoopVectorization]]
-deps = ["DocStringExtensions", "LinearAlgebra", "OffsetArrays", "SIMDPirates", "SLEEFPirates", "UnPack", "VectorizationBase"]
-git-tree-sha1 = "8c19b0c629dae1cea52a990d6303d03fe4328795"
-uuid = "bdcacae8-1622-11e9-2a5c-532679323890"
-version = "0.8.19"
-
-[[METIS_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "3f52ed323683398498ef163a45ce998f1ceca363"
-uuid = "d00139f3-1899-568f-a2f0-47f597d42d70"
-version = "5.1.0+4"
-
-[[MKL_jll]]
-deps = ["IntelOpenMP_jll", "Libdl", "Pkg"]
-git-tree-sha1 = "0ce9a7fa68c70cf83c49d05d2c04d91b47404b08"
-uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
-version = "2020.1.216+0"
-
-[[MLStyle]]
-git-tree-sha1 = "67f9a88611bc79f992aa705d9bbc833a2547dec7"
-uuid = "d8e11817-5142-5d16-987a-aa16d5891078"
-version = "0.3.1"
-
-[[MacroTools]]
-deps = ["Markdown", "Random"]
-git-tree-sha1 = "f7d2e3f654af75f01ec49be82c231c382214223a"
-uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
-version = "0.5.5"
-
-[[Markdown]]
-deps = ["Base64"]
-uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
-
-[[MbedTLS]]
-deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
-git-tree-sha1 = "426a6978b03a97ceb7ead77775a1da066343ec6e"
-uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
-version = "1.0.2"
-
-[[MbedTLS_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "a0cb0d489819fa7ea5f9fa84c7e7eba19d8073af"
-uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
-version = "2.16.6+1"
-
-[[Measures]]
-git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f"
-uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
-version = "0.3.1"
-
-[[Missings]]
-deps = ["DataAPI"]
-git-tree-sha1 = "de0a5ce9e5289f27df672ffabef4d1e5861247d5"
-uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
-version = "0.4.3"
-
-[[Mmap]]
-uuid = "a63ad114-7e13-5084-954f-fe012c677804"
-
-[[ModelingToolkit]]
-deps = ["ArrayInterface", "DiffEqBase", "DiffRules", "Distributed", "DocStringExtensions", "GeneralizedGenerated", "Latexify", "LinearAlgebra", "MacroTools", "NaNMath", "SafeTestsets", "SparseArrays", "SpecialFunctions", "StaticArrays", "TreeViews", "UnPack", "Unitful"]
-git-tree-sha1 = "b3a75a48ddb401d30b2548e72609c27675bde25f"
-uuid = "961ee093-0014-501f-94e3-6117800e7a78"
-version = "3.1.1"
-
-[[MuladdMacro]]
-git-tree-sha1 = "c6190f9a7fc5d9d5915ab29f2134421b12d24a68"
-uuid = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
-version = "0.2.2"
-
-[[MultiScaleArrays]]
-deps = ["DiffEqBase", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "OrdinaryDiffEq", "Random", "RecursiveArrayTools", "SparseDiffTools", "Statistics", "StochasticDiffEq", "TreeViews"]
-git-tree-sha1 = "258f3be6770fe77be8870727ba9803e236c685b8"
-uuid = "f9640e96-87f6-5992-9c3b-0743c6a49ffa"
-version = "1.8.1"
-
-[[MultivariateStats]]
-deps = ["Arpack", "LinearAlgebra", "SparseArrays", "Statistics", "StatsBase"]
-git-tree-sha1 = "352fae519b447bf52e6de627b89f448bcd469e4e"
-uuid = "6f286f6a-111f-5878-ab1e-185364afe411"
-version = "0.7.0"
-
-[[NLSolversBase]]
-deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"]
-git-tree-sha1 = "7c4e66c47848562003250f28b579c584e55becc0"
-uuid = "d41bc354-129a-5804-8e4c-c37616107c6c"
-version = "7.6.1"
-
-[[NLsolve]]
-deps = ["Distances", "LineSearches", "LinearAlgebra", "NLSolversBase", "Printf", "Reexport"]
-git-tree-sha1 = "ea172c86745810136d744fc67650d2e2de669c4f"
-uuid = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
-version = "4.4.0"
-
-[[NNlib]]
-deps = ["BinaryProvider", "Libdl", "LinearAlgebra", "Requires", "Statistics"]
-git-tree-sha1 = "d9f196d911f55aeaff11b11f681b135980783824"
-uuid = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
-version = "0.6.6"
-
-[[NaNMath]]
-git-tree-sha1 = "c84c576296d0e2fbb3fc134d3e09086b3ea617cd"
-uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
-version = "0.3.4"
-
-[[NameResolution]]
-deps = ["DataStructures", "PrettyPrint"]
-git-tree-sha1 = "f4119274d5a410c64a1d9f546312bb6ae54d41c0"
-uuid = "71a1bf82-56d0-4bbc-8a3c-48b961074391"
-version = "0.1.3"
-
-[[NearestNeighbors]]
-deps = ["Distances", "StaticArrays"]
-git-tree-sha1 = "93107e3cdada73d63245ed8170dcae680f0c8fd8"
-uuid = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
-version = "0.4.6"
-
-[[Observables]]
-git-tree-sha1 = "11832878355305984235a2e90d0e3737383c634c"
-uuid = "510215fc-4207-5dde-b226-833fc4488ee2"
-version = "0.3.1"
-
-[[OffsetArrays]]
-git-tree-sha1 = "4ba4cd84c88df8340da1c3e2d8dcb9d18dd1b53b"
-uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "1.1.1"
-
-[[Ogg_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "59cf7a95bf5ac39feac80b796e0f39f9d69dc887"
-uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
-version = "1.3.4+0"
-
-[[OpenBLAS_jll]]
-deps = ["CompilerSupportLibraries_jll", "Libdl", "Pkg"]
-git-tree-sha1 = "0c922fd9634e358622e333fc58de61f05a048492"
-uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.9+5"
-
-[[OpenSSL_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "7aaaded15bf393b5f34c2aad5b765c18d26cb495"
-uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
-version = "1.1.1+4"
-
-[[Optim]]
-deps = ["Compat", "FillArrays", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"]
-git-tree-sha1 = "62054d469d3631960e3f472ceb8624be5b11c34d"
-uuid = "429524aa-4258-5aef-a3af-852621145aeb"
-version = "0.20.6"
-
-[[Opus_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "002c18f222a542907e16c83c64a1338992da7e2c"
-uuid = "91d4177d-7536-5919-b921-800302f37372"
-version = "1.3.1+1"
-
-[[OrderedCollections]]
-git-tree-sha1 = "293b70ac1780f9584c89268a6e2a560d938a7065"
-uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
-version = "1.3.0"
-
-[[OrdinaryDiffEq]]
-deps = ["Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "ExponentialUtilities", "FiniteDiff", "ForwardDiff", "GenericSVD", "LinearAlgebra", "Logging", "MacroTools", "MuladdMacro", "NLsolve", "RecursiveArrayTools", "Reexport", "SparseArrays", "SparseDiffTools", "StaticArrays", "UnPack"]
-git-tree-sha1 = "89f71b92de3ab39d70bd21a2b11346db5aeb9994"
-uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
-version = "5.41.0"
-
-[[PDMats]]
-deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse", "Test"]
-git-tree-sha1 = "b3405086eb6a974eba1958923d46bc0e1c2d2d63"
-uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
-version = "0.10.0"
-
-[[ParameterizedFunctions]]
-deps = ["DataStructures", "DiffEqBase", "Latexify", "LinearAlgebra", "ModelingToolkit", "Reexport"]
-git-tree-sha1 = "458f179b2bcb8071a2dbcf92cafd2c1cbd44e8e1"
-uuid = "65888b18-ceab-5e60-b2b9-181511a3b968"
-version = "5.4.0"
-
-[[Parameters]]
-deps = ["OrderedCollections", "UnPack"]
-git-tree-sha1 = "38b2e970043613c187bd56a995fe2e551821eb4a"
-uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
-version = "0.12.1"
-
-[[Parsers]]
-deps = ["Dates", "Test"]
-git-tree-sha1 = "10134f2ee0b1978ae7752c41306e131a684e1f06"
-uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "1.0.7"
-
-[[Pkg]]
-deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Test", "UUIDs"]
-uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
-
-[[PlotThemes]]
-deps = ["PlotUtils", "Requires", "Statistics"]
-git-tree-sha1 = "c6f5ea535551b3b16835134697f0c65d06c94b91"
-uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
-version = "2.0.0"
-
-[[PlotUtils]]
-deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "Statistics"]
-git-tree-sha1 = "e18e0e51ff07bf92bb7e06dcb9c082a4e125e20c"
-uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
-version = "1.0.5"
-
-[[Plots]]
-deps = ["Base64", "Contour", "Dates", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "GeometryTypes", "JSON", "LinearAlgebra", "Measures", "NaNMath", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs"]
-git-tree-sha1 = "ba747739872a67bc1a8078aec3313bde075b3fb0"
-uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-version = "1.5.5"
-
-[[PoissonRandom]]
-deps = ["Random", "Statistics", "Test"]
-git-tree-sha1 = "44d018211a56626288b5d3f8c6497d28c26dc850"
-uuid = "e409e4f3-bfea-5376-8464-e040bb5c01ab"
-version = "0.4.0"
-
-[[PooledArrays]]
-deps = ["DataAPI"]
-git-tree-sha1 = "b1333d4eced1826e15adbdf01a4ecaccca9d353c"
-uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
-version = "0.5.3"
-
-[[PositiveFactorizations]]
-deps = ["LinearAlgebra", "Test"]
-git-tree-sha1 = "127c47b91990c101ee3752291c4f45640eeb03d1"
-uuid = "85a6dd25-e78a-55b7-8502-1745935b8125"
-version = "0.2.3"
-
-[[PrettyPrint]]
-git-tree-sha1 = "21e1f88cb73589ec39a1181f77eb444b5276151e"
-uuid = "8162dcfd-2161-5ef2-ae6c-7681170c5f98"
-version = "0.1.0"
-
-[[Printf]]
-deps = ["Unicode"]
-uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-
-[[ProgressLogging]]
-deps = ["Logging", "SHA", "UUIDs"]
-git-tree-sha1 = "59398022b661b6fd569f25de6b18fde39843196a"
-uuid = "33c8b6b6-d38a-422a-b730-caa89a2f386c"
-version = "0.1.3"
-
-[[ProgressMeter]]
-deps = ["Distributed", "Printf"]
-git-tree-sha1 = "2de4cddc0ceeddafb6b143b5b6cd9c659b64507c"
-uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
-version = "1.3.2"
-
-[[PyCall]]
-deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"]
-git-tree-sha1 = "3a3fdb9000d35958c9ba2323ca7c4958901f115d"
-uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
-version = "1.91.4"
-
-[[PyPlot]]
-deps = ["Colors", "LaTeXStrings", "PyCall", "Sockets", "Test", "VersionParsing"]
-git-tree-sha1 = "67dde2482fe1a72ef62ed93f8c239f947638e5a2"
-uuid = "d330b81b-6aea-500a-939a-2ce795aea3ee"
-version = "2.9.0"
-
-[[QuadGK]]
-deps = ["DataStructures", "LinearAlgebra"]
-git-tree-sha1 = "0ab8a09d4478ebeb99a706ecbf8634a65077ccdc"
-uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
-version = "2.4.0"
-
-[[REPL]]
-deps = ["InteractiveUtils", "Markdown", "Sockets"]
-uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
-
-[[Random]]
-deps = ["Serialization"]
-uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-
-[[RandomNumbers]]
-deps = ["Random", "Requires"]
-git-tree-sha1 = "441e6fc35597524ada7f85e13df1f4e10137d16f"
-uuid = "e6cf234a-135c-5ec9-84dd-332b85af5143"
-version = "1.4.0"
-
-[[Ratios]]
-git-tree-sha1 = "37d210f612d70f3f7d57d488cb3b6eff56ad4e41"
-uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439"
-version = "0.4.0"
-
-[[RecipesBase]]
-git-tree-sha1 = "58de8f7e33b7fda6ee39eff65169cd1e19d0c107"
-uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
-version = "1.0.2"
-
-[[RecipesPipeline]]
-deps = ["Dates", "PlotUtils", "RecipesBase"]
-git-tree-sha1 = "d2a58b8291d1c0abae6a91489973f8a92bf5c04a"
-uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c"
-version = "0.1.11"
-
-[[RecursiveArrayTools]]
-deps = ["ArrayInterface", "LinearAlgebra", "RecipesBase", "Requires", "StaticArrays", "Statistics", "ZygoteRules"]
-git-tree-sha1 = "0ffe36b65f0fc4967a42a673c1a9ffa65724dee6"
-uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
-version = "2.5.0"
-
-[[RecursiveFactorization]]
-deps = ["LinearAlgebra", "LoopVectorization", "VectorizationBase"]
-git-tree-sha1 = "4ca0bdad1d69abbd59c35af89a9a2ab6cd5ef0f1"
-uuid = "f2c3362d-daeb-58d1-803e-2bc74f2840b4"
-version = "0.1.4"
-
-[[Reexport]]
-deps = ["Pkg"]
-git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0"
-uuid = "189a3867-3050-52da-a836-e630ba90ab69"
-version = "0.2.0"
-
-[[Requires]]
-deps = ["UUIDs"]
-git-tree-sha1 = "d37400976e98018ee840e0ca4f9d20baa231dc6b"
-uuid = "ae029012-a4dd-5104-9daa-d747884805df"
-version = "1.0.1"
-
-[[ResettableStacks]]
-deps = ["StaticArrays"]
-git-tree-sha1 = "d19e9c93de6020a96dbb2820567c78d0ab8f7248"
-uuid = "ae5879a3-cd67-5da8-be7f-38c6eb64a37b"
-version = "1.0.0"
-
-[[Rmath]]
-deps = ["Random", "Rmath_jll"]
-git-tree-sha1 = "86c5647b565873641538d8f812c04e4c9dbeb370"
-uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
-version = "0.6.1"
-
-[[Rmath_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "d76185aa1f421306dec73c057aa384bad74188f0"
-uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
-version = "0.2.2+1"
-
-[[Roots]]
-deps = ["Printf"]
-git-tree-sha1 = "7bd9d7eee602f0413f24c394386a18cb0d515f36"
-uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
-version = "1.0.3"
-
-[[SHA]]
-uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
-
-[[SIMDPirates]]
-deps = ["VectorizationBase"]
-git-tree-sha1 = "ffa6df58c5d7b148e2da027956089ba232f7c81a"
-uuid = "21efa798-c60a-11e8-04d3-e1a92915a26a"
-version = "0.8.17"
-
-[[SLEEFPirates]]
-deps = ["Libdl", "SIMDPirates", "VectorizationBase"]
-git-tree-sha1 = "67ae90a18aa8c22bf159318300e765fbd89ddf6e"
-uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa"
-version = "0.5.5"
-
-[[SafeTestsets]]
-deps = ["Test"]
-git-tree-sha1 = "36ebc5622c82eb9324005cc75e7e2cc51181d181"
-uuid = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
-version = "0.0.1"
-
-[[ScikitLearn]]
-deps = ["Compat", "Conda", "DataFrames", "Distributed", "IterTools", "LinearAlgebra", "MacroTools", "Parameters", "Printf", "PyCall", "Random", "ScikitLearnBase", "SparseArrays", "StatsBase", "VersionParsing"]
-git-tree-sha1 = "b2dbb141575879beb3ad771fb0314a22617586d3"
-uuid = "3646fa90-6ef7-5e7e-9f22-8aca16db6324"
-version = "0.6.2"
-
-[[ScikitLearnBase]]
-deps = ["LinearAlgebra", "Random", "Statistics"]
-git-tree-sha1 = "7877e55c1523a4b336b433da39c8e8c08d2f221f"
-uuid = "6e75b9c4-186b-50bd-896f-2d2496a4843e"
-version = "0.5.0"
-
-[[Serialization]]
-uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
-
-[[SharedArrays]]
-deps = ["Distributed", "Mmap", "Random", "Serialization"]
-uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
-
-[[Showoff]]
-deps = ["Dates"]
-git-tree-sha1 = "e032c9df551fb23c9f98ae1064de074111b7bc39"
-uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
-version = "0.3.1"
-
-[[SimpleTraits]]
-deps = ["InteractiveUtils", "MacroTools"]
-git-tree-sha1 = "2ee666b24ab8be6a922f9d6c11a86e1a703a7dda"
-uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
-version = "0.9.2"
-
-[[Sockets]]
-uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
-
-[[SortingAlgorithms]]
-deps = ["DataStructures", "Random", "Test"]
-git-tree-sha1 = "03f5898c9959f8115e30bc7226ada7d0df554ddd"
-uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
-version = "0.3.1"
-
-[[SparseArrays]]
-deps = ["LinearAlgebra", "Random"]
-uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-
-[[SparseDiffTools]]
-deps = ["Adapt", "ArrayInterface", "Compat", "DataStructures", "FiniteDiff", "ForwardDiff", "LightGraphs", "LinearAlgebra", "Requires", "SparseArrays", "VertexSafeGraphs"]
-git-tree-sha1 = "93666e93899d995ec37abddde4811f533e49c074"
-uuid = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
-version = "1.9.1"
-
-[[SpecialFunctions]]
-deps = ["BinDeps", "BinaryProvider", "Libdl"]
-git-tree-sha1 = "3bdd374b6fd78faf0119b8c5d538788dbf910c6e"
-uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "0.8.0"
-
-[[StaticArrays]]
-deps = ["LinearAlgebra", "Random", "Statistics"]
-git-tree-sha1 = "016d1e1a00fabc556473b07161da3d39726ded35"
-uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "0.12.4"
-
-[[Statistics]]
-deps = ["LinearAlgebra", "SparseArrays"]
-uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-
-[[StatsBase]]
-deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics"]
-git-tree-sha1 = "a6102b1f364befdb05746f386b67c6b7e3262c45"
-uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
-version = "0.33.0"
-
-[[StatsFuns]]
-deps = ["Rmath", "SpecialFunctions"]
-git-tree-sha1 = "04a5a8e6ab87966b43f247920eab053fd5fdc925"
-uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
-version = "0.9.5"
-
-[[StatsPlots]]
-deps = ["Clustering", "DataStructures", "DataValues", "Distributions", "Interpolations", "KernelDensity", "MultivariateStats", "Observables", "Plots", "RecipesBase", "RecipesPipeline", "Reexport", "StatsBase", "TableOperations", "Tables", "Widgets"]
-git-tree-sha1 = "b9b7fff81f573465fcac4685df1497d968537a9e"
-uuid = "f3b207a7-027a-5e70-b257-86293d7955fd"
-version = "0.14.6"
-
-[[SteadyStateDiffEq]]
-deps = ["DiffEqBase", "DiffEqCallbacks", "LinearAlgebra", "NLsolve", "Reexport"]
-git-tree-sha1 = "75f258513b7ef8b235876f4cf146577ffd545094"
-uuid = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
-version = "1.5.1"
-
-[[StochasticDiffEq]]
-deps = ["ArrayInterface", "DataStructures", "DiffEqBase", "DiffEqJump", "DiffEqNoiseProcess", "FillArrays", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "Logging", "MuladdMacro", "NLsolve", "OrdinaryDiffEq", "Random", "RandomNumbers", "RecursiveArrayTools", "Reexport", "SparseArrays", "SparseDiffTools", "StaticArrays", "UnPack"]
-git-tree-sha1 = "a6bd84689c76ad649fcb9c994a7bdb20802846ac"
-uuid = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
-version = "6.24.0"
-
-[[StructArrays]]
-deps = ["Adapt", "DataAPI", "Tables"]
-git-tree-sha1 = "8099ed9fb90b6e754d6ba8c6ed8670f010eadca0"
-uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
-version = "0.4.4"
-
-[[SuiteSparse]]
-deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
-uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
-
-[[SuiteSparse_jll]]
-deps = ["Libdl", "METIS_jll", "OpenBLAS_jll", "Pkg"]
-git-tree-sha1 = "dad9c6f67fa143b7b1935ef3c6231d88e6dcdec2"
-uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c"
-version = "5.4.0+8"
-
-[[Sundials]]
-deps = ["CEnum", "DataStructures", "DiffEqBase", "Libdl", "LinearAlgebra", "Logging", "Reexport", "SparseArrays", "Sundials_jll"]
-git-tree-sha1 = "17fb16310572884dbe2c64c0ce1b1c3ca127662f"
-uuid = "c3572dad-4567-51f8-b174-8c6c989267f4"
-version = "4.2.5"
-
-[[Sundials_jll]]
-deps = ["Libdl", "OpenBLAS_jll", "Pkg", "SuiteSparse_jll"]
-git-tree-sha1 = "4a4ae2ebefa34779552dbe87683c70c856624ec8"
-uuid = "fb77eaff-e24c-56d4-86b1-d163f2edb164"
-version = "5.2.0+0"
-
-[[TableOperations]]
-deps = ["Tables", "Test"]
-git-tree-sha1 = "208630a14884abd110a8f8008b0882f0d0f5632c"
-uuid = "ab02a1b2-a7df-11e8-156e-fb1833f50b87"
-version = "0.2.1"
-
-[[TableTraits]]
-deps = ["IteratorInterfaceExtensions"]
-git-tree-sha1 = "b1ad568ba658d8cbb3b892ed5380a6f3e781a81e"
-uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
-version = "1.0.0"
-
-[[Tables]]
-deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
-git-tree-sha1 = "c45dcc27331febabc20d86cb3974ef095257dcf3"
-uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
-version = "1.0.4"
-
-[[TaylorSeries]]
-deps = ["InteractiveUtils", "LinearAlgebra", "Markdown", "Requires", "SparseArrays"]
-git-tree-sha1 = "192cf2208fd4ffed35d675513c0d5807cfbcaae1"
-uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
-version = "0.10.6"
-
-[[TerminalLoggers]]
-deps = ["LeftChildRightSiblingTrees", "Logging", "Markdown", "Printf", "ProgressLogging", "UUIDs"]
-git-tree-sha1 = "cbea752b5eef52a3e1188fb31580c3e4fa0cbc35"
-uuid = "5d786b92-1e48-4d6f-9151-6b4477ca9bed"
-version = "0.1.2"
-
-[[Test]]
-deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
-uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-
-[[TreeViews]]
-deps = ["Test"]
-git-tree-sha1 = "8d0d7a3fe2f30d6a7f833a5f19f7c7a5b396eae6"
-uuid = "a2a6695c-b41b-5b7d-aed9-dbfdeacea5d7"
-version = "0.3.0"
-
-[[URIParser]]
-deps = ["Unicode"]
-git-tree-sha1 = "53a9f49546b8d2dd2e688d216421d050c9a31d0d"
-uuid = "30578b45-9adc-5946-b283-645ec420af67"
-version = "0.4.1"
-
-[[UUIDs]]
-deps = ["Random", "SHA"]
-uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
-
-[[UnPack]]
-git-tree-sha1 = "e0cb9715adda456f7657e45377fd3063bf87179a"
-uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
-version = "0.1.0"
-
-[[Unicode]]
-uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
-
-[[Unitful]]
-deps = ["ConstructionBase", "LinearAlgebra", "Random"]
-git-tree-sha1 = "a061dada333813818aa7454f93c63a5cab6ea981"
-uuid = "1986cc42-f94f-5a68-af5c-568840ba703d"
-version = "1.3.0"
-
-[[VectorizationBase]]
-deps = ["CpuId", "Libdl", "LinearAlgebra"]
-git-tree-sha1 = "a2a9b59beffc87e593edbefbb55961c39a2dacef"
-uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
-version = "0.12.25"
-
-[[VersionParsing]]
-git-tree-sha1 = "80229be1f670524750d905f8fc8148e5a8c4537f"
-uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
-version = "1.2.0"
-
-[[VertexSafeGraphs]]
-deps = ["LightGraphs"]
-git-tree-sha1 = "b9b450c99a3ca1cc1c6836f560d8d887bcbe356e"
-uuid = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f"
-version = "0.1.2"
-
-[[Widgets]]
-deps = ["Colors", "Dates", "Observables", "OrderedCollections"]
-git-tree-sha1 = "fc0feda91b3fef7fe6948ee09bb628f882b49ca4"
-uuid = "cc8bc4a8-27d6-5769-a93b-9d913e69aa62"
-version = "0.6.2"
-
-[[WoodburyMatrices]]
-deps = ["LinearAlgebra", "SparseArrays"]
-git-tree-sha1 = "28ffe06d28b1ba8fdb2f36ec7bb079fac81bac0d"
-uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6"
-version = "0.5.2"
-
-[[Zlib_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "622d8b6dc0c7e8029f17127703de9819134d1b71"
-uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
-version = "1.2.11+14"
-
-[[Zygote]]
-deps = ["DiffRules", "FFTW", "FillArrays", "ForwardDiff", "IRTools", "InteractiveUtils", "LinearAlgebra", "MacroTools", "NNlib", "NaNMath", "Random", "Requires", "SpecialFunctions", "Statistics", "ZygoteRules"]
-git-tree-sha1 = "f8329b595c465caf3ca87c4f744e6041a4983e43"
-uuid = "e88e6eb3-aa80-5325-afca-941959d7151f"
-version = "0.4.8"
-
-[[ZygoteRules]]
-deps = ["MacroTools"]
-git-tree-sha1 = "b3b4882cc9accf6731a08cc39543fbc6b669dca8"
-uuid = "700de1a5-db45-46bc-99cf-38207098b444"
-version = "0.2.0"
-
-[[libass_jll]]
-deps = ["Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "Libdl", "Pkg", "Zlib_jll"]
-git-tree-sha1 = "027a304b2a90de84f690949a21f94e5ae0f92c73"
-uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"
-version = "0.14.0+2"
-
-[[libfdk_aac_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "480c7ed04f68ea3edd4c757f5db5b6a0a4e0bd99"
-uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280"
-version = "0.1.6+2"
-
-[[libvorbis_jll]]
-deps = ["Libdl", "Ogg_jll", "Pkg"]
-git-tree-sha1 = "6a66f65b5275dfa799036c8a3a26616a0a271c4a"
-uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a"
-version = "1.3.6+4"
-
-[[x264_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "d89346fe63a6465a9f44e958ac0e3d366af90b74"
-uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a"
-version = "2019.5.25+2"
-
-[[x265_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "61324ad346b00a6e541896b94201c9426591e43a"
-uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76"
-version = "3.0.0+1"
diff --git a/examples/Cloudy/Project.toml b/examples/Cloudy/Project.toml
index 53bf842e6..affff83e2 100644
--- a/examples/Cloudy/Project.toml
+++ b/examples/Cloudy/Project.toml
@@ -3,6 +3,7 @@ CalibrateEmulateSample = "95e48a1f-0bec-4818-9538-3db4340308e3"
 Cloudy = "9e3b23bb-e7cc-4b94-886c-65de2234ba87"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 GaussianProcesses = "891a1506-143c-57d2-908e-e1f8e92e6de9"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
@@ -10,3 +11,4 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
+Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
diff --git a/examples/Lorenz/GModel.jl b/examples/Lorenz/GModel.jl
new file mode 100644
index 000000000..e4dcb703c
--- /dev/null
+++ b/examples/Lorenz/GModel.jl
@@ -0,0 +1,423 @@
+module GModel
+
+using DocStringExtensions
+
+using Random
+using Distributions
+using LinearAlgebra
+using DifferentialEquations
+using FFTW
+
+export run_G
+export run_G_ensemble
+export lorenz_forward
+
+"""
+    GSettings{FT<:AbstractFloat, KT, D}
+
+Structure to hold all information to run the forward model G
+
+# Fields
+$(DocStringExtensions.FIELDS)
+"""
+struct GSettings{FT<:AbstractFloat, KT, D}
+    "model output"
+    out::Array{FT, 1}
+    "time period over which to run the model, e.g., `(0, 1)`"
+    tspan::Tuple{FT, FT}
+    "x domain"
+    x::Array{FT, 1}
+    "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
+end
+
+struct LParams
+	# 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}
+
+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
+
+"""
+function run_G_ensemble(params::Array{FT, 2},
+                        settings::LSettings,
+                        rng_seed=42) where {FT<:AbstractFloat}
+    # Initilize ensemble
+    N_ens = size(params, 1) # params is N_ens x N_params
+    if settings.stats_type == 1
+	    nd = 2
+    elseif settings.stats_type == 2
+	    nd = 1 + (2*settings.kmax)
+    elseif settings.stats_type == 3
+	    nd = 3
+    elseif settings.stats_type == 4
+	    nd = Int64(2*(settings.T/settings.Ts/settings.Tfit))
+    elseif settings.stats_type == 5
+	    nd = Int64(settings.T/settings.Ts/settings.Tfit)
+    end
+    g_ens = zeros(N_ens, nd)
+    # Lorenz parameters
+    #Fp = rand(Normal(0.0, 0.01), N); # Initial perturbation
+    F = params[:, 1] # Forcing
+    if settings.dynamics==2
+        A = params[:, 2] # Transience amplitude
+    else
+        A = F .* 0. # Set to zero, it is not used
+    end
+    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) 
+    end
+
+    return g_ens
+end
+
+# Forward pass of the Lorenz 96 model
+# Inputs: settings: structure with stats_type, t_start, T
+# 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
+end
+
+
+# Solve the Lorenz 96 system 
+# Inputs: 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_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
+
+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
+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
+end
+
+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
+end
+
+###############################
+## Extract statistics
+###############################
+function stats(settings, xn, t) 
+# Define averaging indices range
+	indices = findall(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. / 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
+
+
+
+
+## Forward pass of the Lorenz 96 model
+## Inputs: settings: structure with stats_type, t_start, T
+## 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, F) 
+#	# run the Lorenz simulation
+#	xn, t = lorenz_solve(settings, F)
+#	# Define averaging indices range
+#	indices = findall(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. / 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...)
+#	else 
+#		ArgumentError("Setting "*string(settings.stats_type)*" not implemented.")	
+#	end
+#	return gt
+#end
+#
+
+end # module
diff --git a/examples/Lorenz/Lorenz_example.jl b/examples/Lorenz/Lorenz_example.jl
new file mode 100644
index 000000000..d207a4dfa
--- /dev/null
+++ b/examples/Lorenz/Lorenz_example.jl
@@ -0,0 +1,382 @@
+# Import modules
+include("GModel.jl") # Contains Lorenz 96 source code
+
+# Import modules
+using Distributions  # probability distributions and associated functions
+using LinearAlgebra
+using StatsPlots
+using GaussianProcesses
+using Plots
+using Random
+using JLD2
+
+# CES 
+using CalibrateEmulateSample.EnsembleKalmanProcesses
+using CalibrateEmulateSample.GaussianProcessEmulator
+using CalibrateEmulateSample.MarkovChainMonteCarlo
+using CalibrateEmulateSample.Observations
+using CalibrateEmulateSample.Utilities
+using CalibrateEmulateSample.ParameterDistributionStorage
+
+
+rng_seed = 4137
+Random.seed!(rng_seed)
+
+# Output figure save directory
+figure_save_directory = "/home/mhowland/Codes/CESPlots/figures/lorenz_96_dynamic/debug/T_90/"
+data_save_directory = "/home/mhowland/Codes/CESPlots/data/lorenz_96_dynamic/debug/T_90/"
+
+# Governing settings
+# Characteristic time scale
+τc = 5. # 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 = 90. # 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 
+
+
+###
+###  Define the (true) parameters
+###
+# Define the parameters that we want to learn
+F_true = 8. # Mean F
+A_true = 2.5 # Transient F amplitude
+ω_true = 2. * π / (360. / τc) # Frequency of the transient F
+
+if dynamics==2
+    params_true = [F_true, A_true]
+    param_names = ["F", "A"]
+else
+    params_true = [F_true]
+    param_names = ["F"]
+end
+n_param = length(param_names)
+params_true = reshape(params_true, (1,n_param))
+
+println(n_param)
+println(params_true)
+
+
+###
+###  Define the parameter priors
+###
+# Lognormal prior or normal prior?
+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)))
+    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+1.0, A_true+1.0]
+	prior_stds = [2.0, 0.5*A_true]
+	F_prior_mean = F_true+1.0
+	A_prior_mean = 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)
+end
+
+
+###
+###  Define the data from which we want to learn the parameters
+###
+data_names = ["y0", "y1"]
+
+
+###
+###  L96 model settings
+###
+
+# Lorenz 96 model parameters
+# Behavior changes depending on the size of F, N
+N = 36
+dt = 1/64.
+# Start of integration
+t_start = 800.
+# Data collection length
+if dynamics==1
+    T = 2.
+else
+    T = 360. / τc
+end
+# Batch length
+Ts = 5. / τc # Nondimensionalize by L96 timescale
+# Integration length
+Tfit = Ts_days / τc
+# Initial perturbation
+Fp = rand(Normal(0.0, 0.01), N);
+kmax = 1
+# Prescribe variance or use a number of forward passes to define true interval variability
+var_prescribe = false
+# Use CES to learn ω?
+ω_fixed = true
+
+# 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_params = GModel.LParams(F_true, ω_true, A_true)
+
+###
+###  Generate (artificial) truth samples
+###  Note: The observables y are related to the parameters θ by:
+###        y = G(θ) + η
+###
+
+# Lorenz forward
+# Input: params: [N_ens, N_params]
+# Output: gt: [N_ens, N_data]
+# 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=1)
+
+# Compute internal variability covariance
+if var_prescribe==true
+    n_samples = 100
+    yt = zeros(n_samples, length(gt))
+    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))
+    end
+else
+    println("Using truth values to compute covariance")
+    n_samples = 20; 
+    yt = zeros(n_samples, length(gt))
+    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)
+    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=1)
+    println(Γy)
+end
+# Construct observation object
+truth = Observations.Obs(yt, Γy, data_names)
+
+###
+###  Calibrate: Ensemble Kalman Inversion
+###
+
+# L96 settings for the forward model in the EKP
+# Here, the forward model for the EKP settings can be set distinctly from the truth runs
+lorenz_settings_G = lorenz_settings; # initialize to truth settings
+
+# EKP parameters
+log_transform(a::AbstractArray) = log.(a)
+exp_transform(a::AbstractArray) = exp.(a)
+
+N_ens = 50 # number of ensemble members
+N_iter = 5 # number of EKI iterations
+# initial parameters: N_ens x N_params
+initial_params = construct_initial_ensemble(priors, N_ens; rng_seed=rng_seed)
+
+ekiobj = EnsembleKalmanProcesses.EnsembleKalmanProcess(initial_params, 
+						     truth.mean, 
+						     truth.obs_noise_cov,
+						     Inversion())
+
+# EKI iterations
+println("EKP inversion error:")
+err = zeros(N_iter) 
+for i in 1:N_iter
+    if log_normal==false
+        params_i = deepcopy(ekiobj.u[end])
+    else
+        params_i = deepcopy(exp_transform(ekiobj.u[end]))
+    end
+    g_ens = GModel.run_G_ensemble(params_i, lorenz_settings_G)
+    EnsembleKalmanProcesses.update_ensemble!(ekiobj, g_ens) 
+    err[i] = mean((params_true - mean(params_i,dims=1)).^2)
+    println("Iteration: "*string(i)*", Error: "*string(err[i]))
+end
+
+# EKI results: Has the ensemble collapsed toward the truth?
+println("True parameters: ")
+println(params_true)
+
+println("\nEKI results:")
+if log_normal==false
+    println(mean(deepcopy(ekiobj.u[end]), dims=1))
+else
+    println(mean(deepcopy(exp_transform(ekiobj.u[end])), dims=1))
+end
+
+###
+###  Emulate: Gaussian Process Regression
+###
+
+gppackage = GaussianProcessEmulator.GPJL()
+pred_type = GaussianProcessEmulator.YType()
+
+# Construct kernel:
+# Sum kernel consisting of Matern 5/2 ARD kernel, a Squared Exponential Iso 
+# kernel and white noise. Note that the kernels take the signal standard 
+# deviations on a log scale as input.
+len1 = 1.0
+kern1 = SE(len1, 1.0)
+len2 = zeros(3)
+kern2 = Mat52Ard(len2, 0.0)
+white = Noise(log(2.0))
+# # construct kernel
+GPkernel =  kern1 + kern2 + white
+# Get training points from the EKP iteration number in the second input term  
+#u_tp, g_tp = Utilities.extract_GP_tp(ekiobj, N_iter)
+u_tp, g_tp = Utilities.get_training_points(ekiobj, N_iter)
+normalized = true
+
+# Kernel defined above
+#gpobj = GPEmulator.GPObj(u_tp, g_tp, gppackage; GPkernel=GPkernel, 
+#			 obs_noise_cov=Sy, normalized=normalized,
+#                         noise_learn=false, prediction_type=pred_type)
+
+# Default kernel
+gpobj = GaussianProcessEmulator.GaussianProcess(u_tp, g_tp, gppackage; 
+			 obs_noise_cov=Γy, normalized=normalized,
+                         noise_learn=false, prediction_type=pred_type)
+
+# Check how well the Gaussian Process regression predicts on the
+# true parameters
+if log_normal==false
+    y_mean, y_var = GaussianProcessEmulator.predict(gpobj, reshape(params_true, 1, :), 
+                                   transform_to_real=true)
+else
+	y_mean, y_var = GaussianProcessEmulator.predict(gpobj, 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)
+println("GP MSE: ")
+println(mean((truth.mean - vec(y_mean)).^2))
+
+###
+###  Sample: Markov Chain Monte Carlo
+###
+
+# initial values
+u0 = vec(mean(u_tp, dims=1))
+println("initial parameters: ", u0)
+
+# MCMC parameters    
+mcmc_alg = "rwm" # random walk Metropolis
+
+# First let's run a short chain to determine a good step size
+burnin = 0
+step = 0.1 # first guess
+max_iter = 5000
+yt_sample = truth.mean
+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, gpobj)
+
+# Now begin the actual MCMC
+println("Begin MCMC - with step size ", new_step)
+u0 = vec(mean(u_tp, dims=1))
+
+# reset parameters 
+burnin = 2000
+max_iter = 100000
+
+mcmc = MarkovChainMonteCarlo.MCMC(yt_sample, Γy, priors, new_step, u0, max_iter, 
+                    mcmc_alg, burnin, svdflag=true)
+MarkovChainMonteCarlo.sample_posterior!(mcmc, gpobj, max_iter)
+
+println("Posterior size")
+println(size(mcmc.posterior))
+
+posterior = MarkovChainMonteCarlo.get_posterior(mcmc)      
+
+#post_mean = mean(posterior, dims=1)
+#post_cov = cov(posterior, dims=1)
+post_mean = get_mean(posterior)
+post_cov = get_cov(posterior)
+println("post_mean")
+println(post_mean)
+println("post_cov")
+println(post_cov)
+println("D util")
+println(det(inv(post_cov)))
+println(" ")
+
+# Plot the posteriors together with the priors and the true parameter values
+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
+for idx in 1:n_params
+    if idx == 1
+        param = "F"
+	xs = collect(0:0.1:15)
+	xbounds = [0, 15]
+	#xbounds = [5, 11]
+    elseif idx == 2
+        param = "A"
+	xs = collect(0:0.1:5.0)
+	xbounds = [0, 5.0]
+    elseif idx == 3
+        param = "ω"
+	xs = collect(ω_true-0.2:0.01:ω_true+0.2)
+	xbounds = [xs[1], xs[end]]
+    else
+        throw("not implemented")
+    end
+
+    if log_normal==true
+        xs = log.(xs)
+	xbounds = log.(xbounds)
+    end
+
+    label = "true " * param
+    histogram(mcmc.posterior[:, idx], 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, 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)
+
+    title!(param)
+    savefig(save_directory*"posterior_"*param*"_T_"*string(T)*"_w_"*string(ω_true)*".png")
+    #StatsPlots.savefig(save_directory*"posterior_"*param*"_T_"*string(T)*"_w_"*string(ω_true)*".png")
+end
+
+# Save data
+@save data_save_directory*"u_tp.jld2" u_tp
+@save data_save_directory*"g_tp.jld2" g_tp
diff --git a/examples/Lorenz/Project.toml b/examples/Lorenz/Project.toml
new file mode 100644
index 000000000..c9e8dba58
--- /dev/null
+++ b/examples/Lorenz/Project.toml
@@ -0,0 +1,13 @@
+[deps]
+CalibrateEmulateSample = "95e48a1f-0bec-4818-9538-3db4340308e3"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+GaussianProcesses = "891a1506-143c-57d2-908e-e1f8e92e6de9"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
+ScikitLearn = "3646fa90-6ef7-5e7e-9f22-8aca16db6324"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
diff --git a/src/MarkovChainMonteCarlo.jl b/src/MarkovChainMonteCarlo.jl
index 91ebb163b..4adcb7c48 100644
--- a/src/MarkovChainMonteCarlo.jl
+++ b/src/MarkovChainMonteCarlo.jl
@@ -125,8 +125,10 @@ end
 
 function get_posterior(mcmc::MCMC)
     #Return a parameter distributions object
-    posterior_samples = Samples(mcmc.posterior[mcmc.burnin+1:end,:]; params_are_columns=false) 
-    parameter_constraints = get_all_constraints(mcmc.prior) #live in same space as prior
+    parameter_slices = batch(mcmc.prior)
+    posterior_samples = [Samples(mcmc.posterior[mcmc.burnin+1:end,slice]; params_are_columns=false) 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
     posterior_distribution = ParameterDistribution(posterior_samples, parameter_constraints, parameter_names)
     return posterior_distribution
diff --git a/src/ParameterDistribution.jl b/src/ParameterDistribution.jl
index 1789d8872..5cc2996ae 100644
--- a/src/ParameterDistribution.jl
+++ b/src/ParameterDistribution.jl
@@ -18,7 +18,7 @@ export get_name, get_distribution, get_total_dimension, get_dimensions, get_all_
 export sample_distribution
 export no_constraint, bounded_below, bounded_above, bounded
 export transform_constrained_to_unconstrained, transform_unconstrained_to_constrained
-export get_logpdf, get_cov, get_var, get_mean
+export get_logpdf, get_cov, get_var, get_mean, batch
     
 ## Objects
 # for the Distribution
@@ -167,7 +167,7 @@ struct ParameterDistribution{PDType <: ParameterDistributionType, CType <: Const
                                                                        CType <: ConstraintType,
                                                                        ST <: AbstractString}
         
-        parameter_distributions = isa(parameter_distributions, PDType) ? [parameter_distributions] : parameter_distributions
+	parameter_distributions = isa(parameter_distributions, PDType) ? [parameter_distributions] : parameter_distributions
         n_parameter_per_dist = [dimension(pd) for pd in parameter_distributions]
         
         constraints = isa(constraints, Union{<:ConstraintType,Array{<:ConstraintType}}) ? [constraints] : constraints #to calc n_constraints_per_dist
@@ -175,7 +175,7 @@ struct ParameterDistribution{PDType <: ParameterDistributionType, CType <: Const
             
         n_constraints_per_dist = [len(c) for c in constraints]
         n_dists = length(parameter_distributions)
-        n_names = length(names)        
+        n_names = length(names)   
         if !(n_parameter_per_dist == n_constraints_per_dist)
             throw(DimensionMismatch("There must be one constraint per parameter in a distribution, use NoConstraint() type if no constraint is required"))
         elseif !(n_dists == n_names)
@@ -229,7 +229,7 @@ end
 returns the (flattened) array of constraints of the parameter distribution
 """
 function get_all_constraints(pd::ParameterDistribution)
-    return pd.constraints
+    return pd.constraints 
 end
 
 """