Repository: rhayes777/PyAutoFit (python)

Issue:
`LogGaussianPrior` raises Exception
The following script uses ` LogGaussianPrior` and produces an exception:

https://github.com/Jammy2211/autofit_workspace_test/blob/main/scripts/searches/DynestyStatic.py

The exception:

```
2024-02-12 09:53:38,721 - DynestyStatic - INFO - The output path of this fit is /mnt/c/Users/Jammy/Code/PyAuto/autofit_workspace_test/output/searches/DynestyStatic/431e743943ea071dcffa55a0db2535b6
2024-02-12 09:53:38,721 - DynestyStatic - INFO - Outputting pre-fit files (e.g. model.info, visualization).
2024-02-12 09:53:38,948 - DynestyStatic - INFO - Starting new Dynesty non-linear search (no previous samples found).
2024-02-12 09:53:38,993 - autofit.non_linear.initializer - INFO - Generating initial samples of model, which are subject to prior limits and other constraints.
/mnt/c/Users/Jammy/Code/PyAuto/PyAutoFit/autofit/messages/transform.py:195: RuntimeWarning: divide by zero encountered in log
  return self.func(x, *self.args)
1122it [00:04, 241.60it/s, +50 | bound: 40 | nc: 1 | ncall: 4861 | eff(%): 24.361 | loglstar:   -inf < -43.295 <    inf | logz: -62.652 +/-  0.593 | dlogz:  0.001 >  0.059]
2024-02-12 09:53:43,694 - DynestyStatic - INFO - Fit Still Running: Updating results after 10000 iterations (see output folder for latest visualization, samples, etc.)
Traceback (most recent call last):
  File "/mnt/c/Users/Jammy/Code/PyAuto/PyAutoFit/autofit/non_linear/search/abstract_search.py", line 1053, in samples_from
    return self.samples_via_internal_from(
  File "/mnt/c/Users/Jammy/Code/PyAuto/PyAutoFit/autofit/non_linear/search/nest/dynesty/search/abstract.py", line 237, in samples_via_internal_from
    log_prior_list = model.log_prior_list_from(parameter_lists=parameter_lists)
  File "/mnt/c/Users/Jammy/Code/PyAuto/PyAutoFit/autofit/mapper/prior_model/abstract.py", line 1578, in log_prior_list_from
    return [
  File "/mnt/c/Users/Jammy/Code/PyAuto/PyAutoFit/autofit/mapper/prior_model/abstract.py", line 1579, in <listcomp>
    sum(self.log_prior_list_from_vector(vector=vector))
  File "/mnt/c/Users/Jammy/Code/PyAuto/PyAutoFit/autofit/mapper/prior_model/abstract.py", line 1077, in log_prior_list_from_vector
    return list(
  File "/mnt/c/Users/Jammy/Code/PyAuto/PyAutoFit/autofit/mapper/prior_model/abstract.py", line 1079, in <lambda>
    lambda prior_tuple, value: prior_tuple.prior.log_prior_from_value(
  File "/mnt/c/Users/Jammy/Code/PyAuto/PyAutoFit/autofit/mapper/prior/abstract.py", line 188, in __getattr__
    return getattr(self.message, item)
AttributeError: 'TransformedMessage' object has no attribute 'log_prior_from_value'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/mnt/c/Users/Jammy/Code/PyAuto/autofit_workspace_test/scripts/searches/DynestyStatic.py", line 89, in <module>
    result = search.fit(model=model, analysis=analysis)
  File "/mnt/c/Users/Jammy/Code/PyAuto/PyAutoFit/autofit/non_linear/search/abstract_search.py", line 565, in fit
    result = self.start_resume_fit(
  File "/mnt/c/Users/Jammy/Code/PyAuto/PyAutoFit/autofit/non_linear/search/abstract_search.py", line 683, in start_resume_fit
    search_internal = self._fit(
  File "/mnt/c/Users/Jammy/Code/PyAuto/PyAutoFit/autofit/non_linear/search/nest/dynesty/search/abstract.py", line 196, in _fit
    self.perform_update(
  File "/mnt/c/Users/Jammy/Code/PyAuto/PyAutoFit/autofit/non_linear/search/abstract_search.py", line 914, in perform_update
    samples = self.samples_from(model=model, search_internal=search_internal)
  File "/mnt/c/Users/Jammy/Code/PyAuto/PyAutoFit/autofit/non_linear/search/abstract_search.py", line 1057, in samples_from
    return self.samples_via_csv_from(model=model)
  File "/mnt/c/Users/Jammy/Code/PyAuto/PyAutoFit/autofit/non_linear/search/abstract_search.py", line 1080, in samples_via_csv_from
    return self.samples_cls.from_csv(
  File "/mnt/c/Users/Jammy/Code/PyAuto/PyAutoFit/autofit/non_linear/samples/samples.py", line 203, in from_csv
    sample_list = paths.load_samples()
  File "/mnt/c/Users/Jammy/Code/PyAuto/PyAutoFit/autofit/non_linear/paths/directory.py", line 216, in load_samples
    return load_from_table(filename=self._samples_file)
  File "/mnt/c/Users/Jammy/Code/PyAuto/PyAutoFit/autofit/non_linear/samples/sample.py", line 292, in load_from_table
    with open(filename, "r+", newline="") as f:
FileNotFoundError: [Errno 2] No such file or directory: '/mnt/c/Users/Jammy/Code/PyAuto/autofit_workspace_test/output/searches/DynestyStatic/431e743943ea071dcffa55a0db2535b6/files/samples.csv'
```

Pull request:
feature/loguniform
Matt's implementation of LogGaussianPrior.log_prior_from_value

should resolve #932


Patch:
diff --git a/autofit/messages/__init__.py b/autofit/messages/__init__.py
index e16cf974d..ac94276b3 100644
--- a/autofit/messages/__init__.py
+++ b/autofit/messages/__init__.py
@@ -7,6 +7,7 @@ from .normal import (
     NaturalNormal, 
     UniformNormalMessage,
     LogNormalMessage,
+    Log10NormalMessage,
     MultiLogitNormalMessage,
     Log10UniformNormalMessage,
 )
diff --git a/autofit/messages/composed_transform.py b/autofit/messages/composed_transform.py
index 6de8d2eaa..a6ed05fb3 100644
--- a/autofit/messages/composed_transform.py
+++ b/autofit/messages/composed_transform.py
@@ -3,7 +3,7 @@ from typing import Tuple, Optional, Union
 
 import numpy as np
 
-from autofit.messages.abstract import MessageInterface
+from autofit.messages.abstract import MessageInterface, AbstractMessage
 from autofit.messages.transform import AbstractDensityTransform
 
 
@@ -94,6 +94,13 @@ class TransformedMessage(MessageInterface):
         self.lower_limit = lower_limit
         self.upper_limit = upper_limit
 
+        x0, x1 = zip(*base_message._support)
+        z0 = self._inverse_transform(np.array(x0))
+        z1 = self._inverse_transform(np.array(x1))
+        self._support = tuple(zip(z0, z1))
+
+    log_normalisation = AbstractMessage.log_normalisation
+
     def from_natural_parameters(self, new_params, **kwargs):
         return self.with_base(
             self.base_message.from_natural_parameters(new_params, **kwargs,)
@@ -193,10 +200,23 @@ class TransformedMessage(MessageInterface):
             x = _transform.inv_transform(x)
         return x
 
-    def transform_det(self, x):
-        for _transform in self.transforms:
-            x = _transform.log_det(x)
-        return x
+    def _transform_det(self, x):
+        logd = 0
+        for _transform in reversed(self.transforms):
+            x, _logd = _transform.transform_det(x)
+            logd += _logd
+        return x, logd
+    
+    def _transform_det_jac(self, x):
+        logd = 0
+        logd_jacs = []
+        for _transform in reversed(self.transforms):
+            x, _logd, _logd_grad, _jac = _transform.transform_det_jac(x)
+            logd += _logd
+            logd_jacs.append((_logd_grad, _jac))
+
+        return x, logd, logd_jacs
+
 
     def invert_natural_parameters(
         self, natural_parameters: np.ndarray,
@@ -246,14 +266,9 @@ class TransformedMessage(MessageInterface):
     @inverse_transform
     def _sample(self, n_samples) -> np.ndarray:
         return self.base_message._sample(n_samples)
-
-    def _factor(self, _, x: Union[np.ndarray, float],) -> np.ndarray:
-        log_det = self.transform_det(x)
-        x = self._transform(x)
-        eta = self.base_message._broadcast_natural_parameters(x)
-        t = self.base_message.to_canonical_form(x)
-        log_base = self.calc_log_base_measure(x) + log_det
-        return self.base_message.natural_logpdf(eta, t, log_base, self.log_partition)
+    
+    def exp_factor(self, x):
+        return np.exp(np.nan_to_num(self.factor(x), nan=-np.inf))
 
     def factor(self, x: Union[float, np.ndarray]) -> Union[np.ndarray, float]:
         """
@@ -269,7 +284,32 @@ class TransformedMessage(MessageInterface):
         -------
         The probability this value is correct
         """
-        return self._factor(self, x)
+        x, logd = self._transform_det(x)
+        return self.base_message.logpdf(x) + logd
+    
+    def factor_gradient(self, x: Union[float, np.ndarray]) -> Tuple[Union[np.ndarray, float],Union[np.ndarray, float]]:
+        """
+        Call the factor. The closer to the mean a given value is the higher
+        the probability returned.
+
+        Parameters
+        ----------
+        x
+            A value in the space of the transformed message.
+
+        Returns
+        -------
+        The probability this value is correct
+        """
+        x, logd, logd_grad, jacs = self._transform_det_jac(x)
+        logp, grad = self.base_message.logpdf_gradient(x)
+        for jac in reversed(jacs):
+            grad = grad * jac 
+        return logp + logd, grad + logd_grad
+    
+
+
+    
 
     @property
     def multivariate(self):
diff --git a/autofit/messages/normal.py b/autofit/messages/normal.py
index bda48cfb5..8f8870461 100644
--- a/autofit/messages/normal.py
+++ b/autofit/messages/normal.py
@@ -295,6 +295,7 @@ UniformNormalMessage = TransformedMessage(NormalMessage(0, 1), phi_transform)
 Log10UniformNormalMessage = TransformedMessage(UniformNormalMessage, log_10_transform)
 
 LogNormalMessage = TransformedMessage(NormalMessage(0, 1), log_transform)
+Log10NormalMessage = TransformedMessage(NormalMessage(0, 1), log_10_transform)
 
 # Support is the simplex
 MultiLogitNormalMessage = TransformedMessage(
diff --git a/autofit/messages/transform.py b/autofit/messages/transform.py
index a7158b50c..34132027f 100644
--- a/autofit/messages/transform.py
+++ b/autofit/messages/transform.py
@@ -243,9 +243,16 @@ log_transform = FunctionTransform(np.log, np.exp, np.reciprocal, func_grad_hess=
 def _log_10_inv(x):
     return 10**x
 
+LOG10 = np.log(10)
+def _log10_recip(x):
+    return np.reciprocal(x) / LOG10
+
+def log10_3(x):
+    ix = _log10_recip(x)
+    return np.log10(x), ix, -ix/x
 
 # TODO: what should func_grad_hess look like? np.log10, ... ?
-log_10_transform = FunctionTransform(np.log10, _log_10_inv, np.reciprocal)
+log_10_transform = FunctionTransform(np.log10, _log_10_inv, _log10_recip, func_grad_hess=log10_3)
 
 
 def sigmoid(x, scale=1, shift=0):


Test patch:
diff --git a/test_autofit/graphical/functionality/test_messages.py b/test_autofit/graphical/functionality/test_messages.py
index 3ca5d0770..9da8dd84e 100644
--- a/test_autofit/graphical/functionality/test_messages.py
+++ b/test_autofit/graphical/functionality/test_messages.py
@@ -46,15 +46,23 @@ def check_numerical_gradient_hessians(message, x=None):
     res = message.logpdf_gradient(x)
     nres = message.numerical_logpdf_gradient(x)
     for i, (x1, x2) in enumerate(zip(res, nres)):
-        if not np.allclose(x1, x2, rtol=1e-3, atol=1e-2):
-            assert False
+        assert np.allclose(x1, x2, rtol=1e-2, atol=1e-2), (i, x1, x2, message)
 
     res = message.logpdf_gradient_hessian(x)
     nres = message.numerical_logpdf_gradient_hessian(x)
     for i, (x1, x2) in enumerate(zip(res, nres)):
-        assert np.allclose(x1, x2, rtol=1e-3, atol=1e-2), (i, x1, x2, message)
+        assert np.allclose(x1, x2, rtol=1e-2, atol=1e-2), (i, x1, x2, message)
 
 
+N = NormalMessage
+UN = UniformNormalMessage
+SUN = TransformedMessage(UN, LinearShiftTransform(shift=0.3, scale=0.8))
+LN = LogNormalMessage
+MLN = MultiLogitNormalMessage
+# test doubly transformed distributions
+WN = TransformedMessage(
+    NormalMessage(0, 1), transform.log_transform, transform.exp_transform,
+)
 def test_message_norm():
     messages = [
         tuple(map(NormalMessage, [0.5, 0.1], [0.2, 0.3])),
@@ -62,6 +70,7 @@ def test_message_norm():
         tuple(map(GammaMessage, [0.5, 1.1], [0.2, 1.3])),
         tuple(map(GammaMessage, [0.5, 1.1, 2], [0.2, 1.3, 1])),
         tuple(map(BetaMessage, [2.0, 3.2, 1.5], [4.1, 2.3, 3])),
+        tuple(map(NormalMessage, [0.5, 0.1], [0.2, 0.3])),
     ]
     for ms in messages:
         check_log_normalisation(ms)
@@ -70,15 +79,23 @@ def test_message_norm():
             check_numerical_gradient_hessians(m)
 
 
-N = NormalMessage
-UN = UniformNormalMessage
-SUN = TransformedMessage(UN, LinearShiftTransform(shift=0.3, scale=0.8))
-LN = LogNormalMessage
-MLN = MultiLogitNormalMessage
-# test doubly transformed distributions
-WN = TransformedMessage(
-    NormalMessage(0, 1), transform.log_transform, transform.exp_transform,
-)
+def check_transform_norm(dist):
+    norm, err = integrate.quad(dist.exp_factor, *dist._support[0])
+    assert norm == pytest.approx(1, abs=err), dist
+
+
+def test_transform_norm():
+    messages = [
+        tuple(map(UN, [0.5, 0.1, -0.5], [0.2, 0.3, 1.3])),
+        tuple(map(SUN, [0.5, 1.1], [0.2, 1.3])),
+        tuple(map(LN, [0.5, 1.1, 2], [0.2, 1.3, 1])),
+        tuple(map(WN, [2.0, 3.2, 1.5], [4.1, 2.3, 3])),
+    ]
+    for ms in messages:
+        # check_log_normalisation(ms) # TODO get this working
+        for m in ms:
+            check_transform_norm(m)
+            check_numerical_gradient_hessians(m)
 
 
 @pytest.mark.parametrize(
diff --git a/test_autofit/graphical/gaussian/test_prior_factor.py b/test_autofit/graphical/gaussian/test_prior_factor.py
index b80df5680..2e973fa2e 100644
--- a/test_autofit/graphical/gaussian/test_prior_factor.py
+++ b/test_autofit/graphical/gaussian/test_prior_factor.py
@@ -20,10 +20,3 @@ def make_prior_factor(prior):
 
 def test_prior_factor(prior_factor):
     assert prior_factor.prior_model.prior_count == 1
-
-
-def test_log_likelihood_function(prior_factor):
-    instance = prior_factor.prior_model.instance_from_prior_medians()
-    assert prior_factor.log_likelihood_function(instance) == pytest.approx(
-        -0.9189385332046727, 1.0e-8
-    )
diff --git a/test_autofit/test_correspondence.py b/test_autofit/test_correspondence.py
index c7eaccbdc..94b44155f 100644
--- a/test_autofit/test_correspondence.py
+++ b/test_autofit/test_correspondence.py
@@ -53,7 +53,7 @@ def test_to_canonical_form(message):
 
 def test_factor(message):
     x = 15
-    assert message.factor(x) == -4.14213901831845
+    assert np.isclose(message.factor(x), - np.log(20), atol=1e-10)
 
 
 def test_value_for(message, x):