search

CUQI-DTU / CUQIpy

https://cuqi-dtu.github.io/CUQIpy/
Apache License 2.0
48 stars 9 forks source link

Fix exception type raised when gradient is not implemented for UserDefinedDistribution #455

Closed amal-ghamdi closed 5 months ago

amal-ghamdi commented 5 months ago

Description

Calling the BayesainProblem MAP when the prior is UserDefinedDistribution with no gradient fails. It is not supposed to fail but uses an approximate gradient.

Example to reproduce

Provide a code example to reproduce the bug.

# Create the Bayesian model
#--------------------------
import cuqi
import numpy as np
gaussian = cuqi.distribution.Gaussian(np.zeros(2)+0.5, 0.3**2)
print(gaussian.sample())
x_test = cuqi.distribution.UserDefinedDistribution(logpdf_func=gaussian.logpdf, sample_func=gaussian.sample, dim=2)
print(x_test.name)
A_test = cuqi.model.Model(forward=lambda x_test: x_test**2, range_geometry=2, domain_geometry=2)
y_test = cuqi.distribution.Gaussian(A_test(x_test), 0.1**2)
x_test_true = x_test.sample(1)
y_test_data = y_test(x_test=x_test_true).sample()
BP_test = cuqi.problem.BayesianProblem(x_test, y_test)
BP_test.set_data(y_test=y_test_data)
#BP_test.posterior.enable_FD()
map = BP_test.MAP()

output

---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
Cell In [4], [line 16](vscode-notebook-cell:?execution_count=4&line=16)
     [14](vscode-notebook-cell:?execution_count=4&line=14) BP_test.set_data(y_test=y_test_data)
     [15](vscode-notebook-cell:?execution_count=4&line=15) #BP_test.posterior.enable_FD()
---> [16](vscode-notebook-cell:?execution_count=4&line=16) map = BP_test.MAP()

File ~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/problem/_problem.py:284, in BayesianProblem.MAP(self, disp, x0)
    [281](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/problem/_problem.py:281)     solver_info = {"solver": "direct"}
    [283](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/problem/_problem.py:283) else: # If no specific implementation exists, use numerical optimization.
--> [284](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/problem/_problem.py:284)     x_MAP, solver_info = self._solve_max_point(self.posterior, disp=disp, x0=x0)
    [286](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/problem/_problem.py:286) # Wrap the result in a CUQIarray and add solver info
    [287](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/problem/_problem.py:287) x_MAP = cuqi.array.CUQIarray(x_MAP, geometry=self.posterior.geometry)

File ~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/problem/_problem.py:651, in BayesianProblem._solve_max_point(self, density, disp, x0)
    [648](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/problem/_problem.py:648)     if disp: print("Optimizing with approximate gradients.") 
    [650](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/problem/_problem.py:650) # Compute point estimate
--> [651](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/problem/_problem.py:651) if self._check_posterior(self, CMRF, must_have_gradient=True): # Use L-BFGS-B for CMRF prior as it has better performance for this multi-modal posterior
    [652](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/problem/_problem.py:652)     if disp: print(f"Using scipy.optimize.L_BFGS_B on negative log of {density.__class__.__name__}")
    [653](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/problem/_problem.py:653)     if disp: print("x0: ones vector")

File ~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/problem/_problem.py:713, in BayesianProblem._check_posterior(posterior, prior_type, likelihood_type, model_type, max_dim, must_have_gradient)
    [711](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/problem/_problem.py:711) if must_have_gradient:
    [712](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/problem/_problem.py:712)     try: 
--> [713](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/problem/_problem.py:713)         posterior.prior.gradient(np.zeros(posterior.prior.dim))
    [714](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/problem/_problem.py:714)         posterior.likelihood.gradient(np.zeros(posterior.likelihood.dim))
    [715](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/problem/_problem.py:715)         G = True

File ~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/density/_density.py:102, in Density.gradient(self, *args, **kwargs)
     [98](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/density/_density.py:98)     return cuqi.utilities.approx_gradient(
     [99](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/density/_density.py:99)         self.logd, *args, **kwargs, epsilon=self.FD_epsilon)
    [101](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/density/_density.py:101) # Otherwise use the implemented gradient
--> [102](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/density/_density.py:102) return self._gradient(*args, **kwargs)

File ~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/distribution/_custom.py:59, in UserDefinedDistribution._gradient(self, x)
     [57](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/distribution/_custom.py:57)     return self.gradient_func(x)
     [58](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/distribution/_custom.py:58) else:
---> [59](https://file+.vscode-resource.vscode-cdn.net/Users/amal/Documents/research_code/CUQI-DTU/Collab-BrainEfflux/code/ear_aqueducts/~/Documents/research_code/CUQI-DTU/CUQIpy/cuqi/distribution/_custom.py:59)     raise Exception("gradient_func is not defined.")

Exception: gradient_func is not defined.