Mean of GMRF seems not properly set #627
Description
Description
As @jeverink reported, for the demo in readme, if we switch LMRF to GMRF, it gives a "dimension mismatch" error. This seems due to that the mean of GMRF is not properly set.
Example to reproduce
Provide a code example to reproduce the bug
# Imports
import matplotlib.pyplot as plt
from cuqi.testproblem import Deconvolution2D
from cuqi.distribution import Gaussian, LMRF, Gamma, GMRF
from cuqi.problem import BayesianProblem
# Step 1: Set up forward model and data, y = Ax
A, y_data, info = Deconvolution2D(dim=256, phantom="cookie").get_components()
# Step 2: Define distributions for parameters
d = Gamma(1, 1e-4)
s = Gamma(1, 1e-4)
x = GMRF(0, lambda d: 1/d, geometry=A.domain_geometry)
y = Gaussian(A@x, lambda s: 1/s)
# Step 3: Combine into Bayesian Problem and sample posterior
BP = BayesianProblem(y, x, d, s)
BP.set_data(y=y_data)
samples = BP.sample_posterior(200)File ~/projects/CUQIpy/cuqi/problem/_problem.py:338, in BayesianProblem.sample_posterior(self, Ns, Nb, callback, experimental)
[335](CUQIpy/~/projects/CUQIpy/cuqi/problem/_problem.py:335) # If target is a joint distribution, try Gibbs sampling
[336](CUQIpy/~/projects/CUQIpy/cuqi/problem/_problem.py:336) # This is still very experimental!
[337](CUQIpy/~/projects/CUQIpy/cuqi/problem/_problem.py:337) if isinstance(self._target, JointDistribution):
--> [338](CUQIpy/~/projects/CUQIpy/cuqi/problem/_problem.py:338) return self._sampleGibbs(Ns, Nb, callback=callback, experimental=experimental)
[340](CUQIpy/~/projects/CUQIpy/cuqi/problem/_problem.py:340) # For Gaussian small-scale we can use direct sampling
[341](CUQIpy/~/projects/CUQIpy/cuqi/problem/_problem.py:341) if self._check_posterior(self, Gaussian, Gaussian, LinearModel, config.MAX_DIM_INV) and not self._check_posterior(self, GMRF):
File ~/projects/CUQIpy/cuqi/problem/_problem.py:888, in BayesianProblem._sampleGibbs(self, Ns, Nb, callback, experimental)
[885](CUQIpy/~/projects/CUQIpy/cuqi/problem/_problem.py:885) sampling_strategy = self._determine_sampling_strategy()
[887](CUQIpy/~/projects/CUQIpy/cuqi/problem/_problem.py:887) sampler = cuqi.sampler.Gibbs(self._target, sampling_strategy)
--> [888](CUQIpy/~/projects/CUQIpy/cuqi/problem/_problem.py:888) samples = sampler.sample(Ns, Nb)
[890](CUQIpy/~/projects/CUQIpy/cuqi/problem/_problem.py:890) # Print timing
[891](CUQIpy/~/projects/CUQIpy/cuqi/problem/_problem.py:891) print('Elapsed time:', time.time() - ti)
File ~/projects/CUQIpy/cuqi/sampler/_gibbs.py:106, in Gibbs.sample(self, Ns, Nb)
[104](CUQIpy/~/projects/CUQIpy/cuqi/sampler/_gibbs.py:104) # Sample tuning phase
[105](CUQIpy/~/projects/CUQIpy/cuqi/sampler/_gibbs.py:105) for i in range(at_Nb, at_Nb+Nb):
...
--> [553](CUQIpy/~/miniconda3/envs/cuqipydev/lib/python3.12/site-packages/scipy/sparse/_base.py:553) raise ValueError('dimension mismatch')
[555](CUQIpy/~/miniconda3/envs/cuqipydev/lib/python3.12/site-packages/scipy/sparse/_base.py:555) result = self._mul_vector(np.ravel(other))
[557](CUQIpy/~/miniconda3/envs/cuqipydev/lib/python3.12/site-packages/scipy/sparse/_base.py:557) if isinstance(other, np.matrix):
ValueError: dimension mismatchand x.mean is set to array([0]) while x.dim is 65536.