CUQI-DTU / CUQIpy-FEniCS

CUQIpy plug-in for FEniCS PDE modeling
GNU General Public License v3.0
4 stars 2 forks source link

bug in sampling with NUTS for specific setup of Gaussian prior #63

Open amal-ghamdi opened 11 months ago

amal-ghamdi commented 11 months ago

The Gaussian prior setup Gaussian(np.zeros(21*21), 1) works while the Gaussian prior setup Gaussian(0, np.ones(21*21)) fails in the scenario below. This is tested against dev #62 and sprint20_UQ_plot in CUQIpy. We should try it with the latest version after these two branches are merged.

This test below failed for the cases 2 and 3,

@pytest.mark.parametrize(
    "exactSolution, noise_level, field_type, field_params, mapping," +
    "prior, case",
    [(None, 0.02, None, None, None, None, "valid"),
     (None, 0.02, None, None, None, None, "unknown_field"),
     (None, 0.05, None, None, "exponential", Gaussian(0, np.ones(21*21)), "valid"),
     (lambda x: x[0]+x[1]+0.1, 0.02, None, None, "exponential",
    Gaussian(0, np.ones(21*21), name='x'), "valid"),
     (None, 0.05, None, None, lambda m: m+1, None, "valid"),
     (None, 0.1, "KL", None, "exponential", None, "valid"),
     (np.random.randn(10), 0.02, "KL",
      {'length_scale': 0.2, 'num_terms': 10}, "exponential", None, "valid")])
def test_FEniCSPoisson2D_setup(exactSolution, noise_level,
                               field_type, field_params, mapping,
                               prior, case):
    """Test creating a FEniCSPoisson2D testproblem with different 
    parametrization, exact solution and noise level"""
    np.random.seed(0)

    if case == "valid":
        testproblem = FEniCSPoisson2D(
            (20, 20),
            exactSolution=exactSolution,
            noise_level=noise_level,
            field_type=field_type,
            field_params=field_params,
            mapping=mapping,
            prior=prior)

        testproblem.data.plot()
        testproblem.exactData.plot()
        testproblem.exactSolution.plot()
        testproblem.UQ(Ns=20, percent=90)
    .
    .
    .

Specifics of test message

../CUQIpy/cuqi/problem/_problem.py:397: in UQ
    samples = self.sample_posterior(Ns)
../CUQIpy/cuqi/problem/_problem.py:332: in sample_posterior
    elif hasattr(self.prior,"sqrtprecTimesMean") and hasattr(self.likelihood.distribution,"sqrtprec") and isinstance(self.model,LinearModel):
../CUQIpy/cuqi/distribution/_gaussian.py:238: in sqrtprecTimesMean
    return (self.sqrtprec@self.mean).flatten()
/Users/amal/opt/miniconda3/envs/fenicsproject/lib/python3.10/site-packages/scipy/sparse/base.py:560: in __matmul__
    return self.__mul__(other)

To make the test pass. I used different representation of the prior: e.g. use Gaussian(np.zeros(21*21), 1) instead of Gaussian(0, np.ones(21*21))

nabriis commented 4 months ago

@amal-ghamdi. Is this still an issue? I can't see an obvious issue with sqrtprecTimesMean