Bugfix parameter name bug for MV distributions, cap mac-os tests to j…

…ulia 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

@@ -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: ""

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,
     )

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...)