From b6a80c9a456ef3df7b70e8c6635e8dcdea920c1e Mon Sep 17 00:00:00 2001 From: Oliver Dunbar <47412152+odunbar@users.noreply.github.com> Date: Mon, 29 Jan 2024 11:16:00 -0800 Subject: [PATCH] Bugfix parameter name bug for MV distributions, cap mac-os tests to julia v1.9.0 (#277) * sorted parameter name multiplicity for MV distributions * add test for MV case, and comment redundant if-else branches * sqrt bug * missed codecov * check Mac test with Julia 1.9.3 * wrong version * try 1.8.5 macos * try 1.9.2 macos * try 1.9.0 * 1.9.1 try * latest working version 1.9.0` --- .github/workflows/Tests.yml | 6 ++-- src/MarkovChainMonteCarlo.jl | 27 ++++++++++++++---- test/MarkovChainMonteCarlo/runtests.jl | 38 ++++++++++++++++++++++++++ 3 files changed, 62 insertions(+), 9 deletions(-) diff --git a/.github/workflows/Tests.yml b/.github/workflows/Tests.yml index 2e6854f91..53e2c3c99 100644 --- a/.github/workflows/Tests.yml +++ b/.github/workflows/Tests.yml @@ -26,7 +26,7 @@ jobs: uses: julia-actions/setup-julia@v1 with: version: 1 - + - name: Install Julia Project Packages # we add this ENV varaible to force PyCall to download and use Conda rather than # the system python (default on Linux), see the PyCall documentation @@ -80,8 +80,8 @@ jobs: - name: Set up Julia uses: julia-actions/setup-julia@v1 with: - version: 1 - + version: 1.9.0 + # later versions causes hanging in package building - name: Install Julia Project Packages env: PYTHON: "" diff --git a/src/MarkovChainMonteCarlo.jl b/src/MarkovChainMonteCarlo.jl index e7a862c28..3db86e50e 100644 --- a/src/MarkovChainMonteCarlo.jl +++ b/src/MarkovChainMonteCarlo.jl @@ -96,8 +96,8 @@ AdvancedMH.logratio_proposal_density( function _get_proposal(prior::ParameterDistribution) # *only* use covariance of prior, not full distribution - Σ = ParameterDistributions.cov(prior) - return AdvancedMH.RandomWalkProposal(MvNormal(zeros(size(Σ)[1]), Σ)) + Σsqrt = sqrt(ParameterDistributions.cov(prior)) # rt_cov * MVN(0,I) avoids the posdef errors for MVN in Julia Distributions + return AdvancedMH.RandomWalkProposal(Σsqrt * MvNormal(zeros(size(Σsqrt)[1]), I)) end """ @@ -302,6 +302,8 @@ function AbstractMCMC.bundle_samples( # 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))] + # elseif length(param_names) < length(keys(ts[1].params))# in case bug with MV names, Chains still needs one name per dist. + # 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) @@ -310,9 +312,9 @@ function AbstractMCMC.bundle_samples( # Bundle everything up and return a MCChains.Chains struct. return MCMCChains.Chains( - vals, - vcat(param_names, internal_names), - (parameters = param_names, internals = internal_names); + vals, # current state information as vec-of-vecs + vcat(param_names, internal_names), # parameter names which get converted to symbols + (parameters = param_names, internals = internal_names); # name map (one needs to be called parameters = ...) start = discard_initial + 1, thin = thinning, ) @@ -350,6 +352,8 @@ function AbstractMCMC.bundle_samples( # Check if we received any parameter names. if ismissing(param_names) param_names = [Symbol(:param_, i) for i in 1:length(keys(ts[1][1].params))] + # elseif length(param_names) < length(keys(ts[1][1].params)) # in case bug with MV names, Chains still needs one name per dist. + # param_names = [Symbol(:param_, i) for i in 1:length(keys(ts[1][1].params))] else # Generate new array to be thread safe. param_names = Symbol.(param_names) @@ -427,9 +431,20 @@ function MCMCWrapper( obs_sample = to_decorrelated(obs_sample, em) log_posterior_map = EmulatorPosteriorModel(prior, em, obs_sample) mh_proposal_sampler = MetropolisHastingsSampler(mcmc_alg, prior) + + # parameter names are needed in every dimension in a MCMCChains object needed for diagnostics + # so create the duplicates here + dd = get_dimensions(prior) + if all(dd .== 1) # i.e if dd == [1, 1, 1, 1, 1], => all params are univariate + param_names = get_name(prior) + else # else use multiplicity to get still informative parameter names + pn = get_name(prior) + param_names = reduce(vcat, [(pn[k] * "_") .* map(x -> string(x), 1:dd[k]) for k in 1:length(pn)]) + end + sample_kwargs = (; # set defaults here :init_params => deepcopy(init_params), - :param_names => get_name(prior), + :param_names => param_names, :discard_initial => burnin, :chain_type => MCMCChains.Chains, ) diff --git a/test/MarkovChainMonteCarlo/runtests.jl b/test/MarkovChainMonteCarlo/runtests.jl index 3b1c1ec75..de8ce8cb5 100644 --- a/test/MarkovChainMonteCarlo/runtests.jl +++ b/test/MarkovChainMonteCarlo/runtests.jl @@ -34,6 +34,33 @@ function test_prior() return ParameterDistribution(prior_dist, prior_constraint, prior_name) end +function test_prior_mv() + ### Define prior + return constrained_gaussian("u_mv", -1.0, 6.0, -Inf, Inf, repeats = 10) +end + +function test_data_mv(; rng_seed = 41, n = 20, var_y = 0.05, input_dim = 10, rest...) + # Seed for pseudo-random number generator + rng = Random.MersenneTwister(rng_seed) + n = 40 # number of training points + x = 1 / input_dim * π * rand(rng, Float64, (input_dim, n)) # predictors/features: 1 × n + σ2_y = reshape([var_y], 1, 1) + y = sin.(norm(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 +end +function test_gp_mv(y, σ2_y, iopairs::PairedDataContainer; norm_factor = nothing) + gppackage = GPJL() + pred_type = YType() + # Construct kernel: + # Squared exponential kernel (note that hyperparameters are on log scale) + # with observational noise + gp = GaussianProcess(gppackage; noise_learn = true, prediction_type = pred_type) + em = Emulator(gp, iopairs; obs_noise_cov = σ2_y) + Emulators.optimize_hyperparameters!(em) + return em +end + function test_gp_1(y, σ2_y, iopairs::PairedDataContainer; norm_factor = nothing) gppackage = GPJL() pred_type = YType() @@ -117,6 +144,17 @@ end @test isapprox(test_obs, (obs_sample ./ sqrt(σ2_y[1, 1])); atol = 1e-2) end + @testset "MV priors" begin + # 10D dist with 1 name, just build the wrapper for test + prior_mv = test_prior_mv() + y_mv, σ2_y_mv, iopairs_mv, rng_mv = test_data(input_dim = 10) + init_params = repeat([0.0], 10) + obs_sample = obs_sample # scalar or Vector -> Vector + em_mv = test_gp_mv(y_mv, σ2_y_mv, iopairs_mv) + mcmc = MCMCWrapper(RWMHSampling(), obs_sample, prior_mv, em_mv; init_params = init_params) + + end + @testset "Sine GP & RW Metropolis" begin em_1 = test_gp_1(y, σ2_y, iopairs) new_step, posterior_mean_1 = mcmc_test_template(prior, σ2_y, em_1; mcmc_params...)