Repository: lenstronomy/lenstronomy

Issue:
Easy check of parameter bounds
It would be nice to get a quick alert if any of the parameters for the PSO are initialized outside of the parameter bounds.

Hints from issue comments:
@sibirrer: I close this issue for now @amn3142. Please reopen if needed.

Solution patch:
diff --git a/lenstronomy/Sampling/Likelihoods/position_likelihood.py b/lenstronomy/Sampling/Likelihoods/position_likelihood.py
index 25471d99..59c23a39 100644
--- a/lenstronomy/Sampling/Likelihoods/position_likelihood.py
+++ b/lenstronomy/Sampling/Likelihoods/position_likelihood.py
@@ -138,8 +138,7 @@ class PositionLikelihood(object):
         return logL
 
     def check_additional_images(self, kwargs_ps, kwargs_lens):
-        """Checks whether additional images have been found and placed in kwargs_ps of
-        the first point source model.
+        """Checks whether additional images have been found and placed in kwargs_ps.
 
         :param kwargs_ps: point source kwargs
         :param kwargs_lens: lens model keyword arguments
@@ -168,6 +167,7 @@ class PositionLikelihood(object):
         :return: log likelihood of the astrometirc correction between predicted image
             positions and model placement of the point sources
         """
+        # TODO: make it compatible with multiple source instances
         if not len(kwargs_ps) > 0:
             return 0
         if "ra_image" not in kwargs_ps[0]:
@@ -267,9 +267,9 @@ class PositionLikelihood(object):
                         if delta[0] ** 2 + delta[1] ** 2 > hard_bound_rms**2:
                             if verbose is True:
                                 print(
-                                    "Image positions do not match to the same source position to the required "
+                                    "Image positions of image %s of model %s do not match to the same source position to the required "
                                     "precision. Achieved: %s, Required: %s."
-                                    % (delta, hard_bound_rms)
+                                    % (i, k, delta, hard_bound_rms)
                                 )
                             logL -= 10**3
                     try:
diff --git a/lenstronomy/Sampling/Likelihoods/time_delay_likelihood.py b/lenstronomy/Sampling/Likelihoods/time_delay_likelihood.py
index 7cc8bc7c..ec8ed7a6 100644
--- a/lenstronomy/Sampling/Likelihoods/time_delay_likelihood.py
+++ b/lenstronomy/Sampling/Likelihoods/time_delay_likelihood.py
@@ -38,10 +38,13 @@ class TimeDelayLikelihood(object):
         self._pointSource = point_source_class
 
     def logL(self, kwargs_lens, kwargs_ps, kwargs_cosmo):
-        """Routine to compute the log likelihood of the time delay distance :param
-        kwargs_lens: lens model kwargs list :param kwargs_ps: point source kwargs list
-        :param kwargs_cosmo: cosmology and other kwargs :return: log likelihood of the
-        model given the time delay data."""
+        """Routine to compute the log likelihood of the time-delay distance.
+
+        :param kwargs_lens: lens model kwargs list
+        :param kwargs_ps: point source kwargs list
+        :param kwargs_cosmo: cosmology and other kwargs
+        :return: log likelihood of the model given the time delay data.
+        """
         x_pos, y_pos = self._pointSource.image_position(
             kwargs_ps=kwargs_ps, kwargs_lens=kwargs_lens, original_position=True
         )
diff --git a/lenstronomy/Sampling/likelihood.py b/lenstronomy/Sampling/likelihood.py
index 8d7c61b2..d4e4acfa 100644
--- a/lenstronomy/Sampling/likelihood.py
+++ b/lenstronomy/Sampling/likelihood.py
@@ -345,7 +345,7 @@ class LikelihoodModule(object):
                 args, self._lower_limit, self._upper_limit, verbose=verbose
             )
             if bound_hit is True:
-                return -(10**15)
+                return -(10**18)
         return self.log_likelihood(kwargs_return, verbose=verbose)
 
     def log_likelihood(self, kwargs_return, verbose=False):
diff --git a/lenstronomy/Sampling/parameters.py b/lenstronomy/Sampling/parameters.py
index de2b82a7..5fbc0489 100644
--- a/lenstronomy/Sampling/parameters.py
+++ b/lenstronomy/Sampling/parameters.py
@@ -985,6 +985,23 @@ class Param(object):
         else:
             return 0
 
+    @property
+    def fixed_kwargs_list(self):
+        """
+
+        :return: list of fixed keyword arguments
+        """
+
+        return (
+            self.lensParams.kwargs_fixed,
+            self.sourceParams.kwargs_fixed,
+            self.lensLightParams.kwargs_fixed,
+            self.pointSourceParams.kwargs_fixed,
+            self.specialParams.kwargs_fixed,
+            self.extinctionParams.kwargs_fixed,
+            self.tracerSourceParams.kwargs_fixed,
+        )
+
     def print_setting(self):
         """Prints the setting of the parameter class.
 
diff --git a/lenstronomy/Sampling/special_param.py b/lenstronomy/Sampling/special_param.py
index 90e5855f..fe33d739 100644
--- a/lenstronomy/Sampling/special_param.py
+++ b/lenstronomy/Sampling/special_param.py
@@ -290,7 +290,7 @@ class SpecialParam(object):
 
         if kwargs_fixed is None:
             kwargs_fixed = {}
-        self._kwargs_fixed = kwargs_fixed
+        self.kwargs_fixed = kwargs_fixed
 
         if kwargs_lower is None:
             kwargs_lower = {}
@@ -317,13 +317,13 @@ class SpecialParam(object):
                 self._param_groups,
                 args,
                 i,
-                kwargs_fixed=self._kwargs_fixed,
+                kwargs_fixed=self.kwargs_fixed,
                 kwargs_lower=self.lower_limit,
                 kwargs_upper=self.upper_limit,
             )
         else:
             result = ModelParamGroup.compose_get_params(
-                self._param_groups, args, i, kwargs_fixed=self._kwargs_fixed
+                self._param_groups, args, i, kwargs_fixed=self.kwargs_fixed
             )
         return result
 
@@ -334,7 +334,7 @@ class SpecialParam(object):
         :return: argument list of the sampled parameters extracted from kwargs_special
         """
         return ModelParamGroup.compose_set_params(
-            self._param_groups, kwargs_special, kwargs_fixed=self._kwargs_fixed
+            self._param_groups, kwargs_special, kwargs_fixed=self.kwargs_fixed
         )
 
     def num_param(self):
@@ -343,7 +343,7 @@ class SpecialParam(object):
         :return: integer, number of free parameters sampled (and managed) by this class, parameter names (list of strings)
         """
         return ModelParamGroup.compose_num_params(
-            self._param_groups, kwargs_fixed=self._kwargs_fixed
+            self._param_groups, kwargs_fixed=self.kwargs_fixed
         )
 
     @property
diff --git a/lenstronomy/SimulationAPI/observation_api.py b/lenstronomy/SimulationAPI/observation_api.py
index d66b235f..0104dd5b 100644
--- a/lenstronomy/SimulationAPI/observation_api.py
+++ b/lenstronomy/SimulationAPI/observation_api.py
@@ -328,10 +328,12 @@ class SingleBand(Instrument, Observation):
         return noise
 
     def estimate_noise(self, image):
-        """
+        """Estimates noise per pixel on the basis that the counts in the pixel come from
+        Poisson statistics of astrophysical objects.
 
-        :param image: noisy data, background subtracted
-        :return: estimated noise map  sqrt(variance) for each pixel as estimated from the instrument and observation
+        :param image: data, background subtracted
+        :return: estimated noise map sqrt(variance) for each pixel as estimated from the
+            instrument and observation
         """
         return np.sqrt(self.background_noise**2 + self.flux_noise(image) ** 2)
 
diff --git a/lenstronomy/Workflow/fitting_sequence.py b/lenstronomy/Workflow/fitting_sequence.py
index 86cf24c5..7f469b18 100644
--- a/lenstronomy/Workflow/fitting_sequence.py
+++ b/lenstronomy/Workflow/fitting_sequence.py
@@ -88,8 +88,23 @@ class FittingSequence(object):
             fitting_type = fitting[0]
             kwargs = fitting[1]
 
+            if fitting_type in [
+                "PSO",
+                "SIMPLEX",
+                "MCMC",
+                "emcee",
+                "zeus",
+                "Cobaya",
+                "dynesty",
+                "dyPolyChord",
+                "MultiNest",
+                "nested_sampling",
+                "Nautilus",
+            ]:
+                self._updateManager.check_initial_state()
             if fitting_type == "restart":
                 self._updateManager.set_init_state()
+                self._updateManager.check_initial_state()
 
             elif fitting_type == "update_settings":
                 self.update_settings(**kwargs)
@@ -193,7 +208,8 @@ class FittingSequence(object):
     def best_fit(self, bijective=False):
         """
 
-        :param bijective: bool, if True, the mapping of image2source_plane and the mass_scaling parameterisation are inverted. If you do not use those options, there is no effect.
+        :param bijective: bool, if True, the mapping of image2source_plane and the mass_scaling parameterisation
+         are inverted. If you do not use those options, there is no effect.
         :return: best fit model of the current state of the FittingSequence class
         """
 
@@ -207,17 +223,19 @@ class FittingSequence(object):
         """
         self._updateManager.update_param_state(**kwargs_update)
 
-    @property
-    def best_fit_likelihood(self):
+    def best_fit_likelihood(self, verbose=False):
         """Returns the log likelihood of the best fit model of the current state of this
         class.
 
+        :param verbose: bool, if True, prints likelihood statements
         :return: log likelihood, float
         """
         kwargs_result = self.best_fit(bijective=True)
         param_class = self.param_class
         likelihoodModule = self.likelihoodModule
-        logL = likelihoodModule.logL(param_class.kwargs2args(**kwargs_result))
+        logL = likelihoodModule.logL(
+            param_class.kwargs2args(**kwargs_result), verbose=verbose
+        )
         return logL
 
     @property
@@ -230,7 +248,9 @@ class FittingSequence(object):
         num_param_nonlinear = self.param_class.num_param()[0]
         num_param_linear = self.param_class.num_param_linear()
         num_param = num_param_nonlinear + num_param_linear
-        bic = analysis_util.bic_model(self.best_fit_likelihood, num_data, num_param)
+        bic = analysis_util.bic_model(
+            self.best_fit_likelihood(verbose=False), num_data, num_param
+        )
         return bic
 
     @property
diff --git a/lenstronomy/Workflow/update_manager.py b/lenstronomy/Workflow/update_manager.py
index 42929c71..58f65db8 100644
--- a/lenstronomy/Workflow/update_manager.py
+++ b/lenstronomy/Workflow/update_manager.py
@@ -1,4 +1,5 @@
 import copy
+import numpy as np
 from lenstronomy.Sampling.parameters import Param
 
 __all__ = ["UpdateManager"]
@@ -171,6 +172,23 @@ class UpdateManager(object):
             "kwargs_tracer_source": self._tracer_source_init,
         }
 
+    @property
+    def _init_kwargs(self):
+        """Keyword arguments for all model components of the initial mean model
+        proposition in the sampling.
+
+        :return: list of keyword arguments
+        """
+        return (
+            self._lens_init,
+            self._source_init,
+            self._lens_light_init,
+            self._ps_init,
+            self._special_init,
+            self._extinction_init,
+            self._tracer_source_init,
+        )
+
     @property
     def sigma_kwargs(self):
         """
@@ -223,6 +241,63 @@ class UpdateManager(object):
             self._tracer_source_fixed,
         )
 
+    @property
+    def fixed_kwargs_updated_list(self):
+        """
+
+        :return: list of fixed arguments including all settings
+        """
+        return self.param_class.fixed_kwargs_list
+
+    def check_initial_state(self):
+        """Checks if initial state is within the bounds in all parameters returns a
+        warning for specific arguments being out of the bounds.
+
+        #TODO: checks whether parameters match the model definitions
+
+        :return:
+        """
+        _model_class_list = [
+            "kwargs_lens",
+            "kwargs_source",
+            "kwargs_lens_light",
+            "kwargs_ps",
+            "kwargs_special",
+            "kwargs_extinction",
+            "kwargs_tracer_source",
+        ]
+        init_kwargs_list = self._init_kwargs
+        lower_kwargs_list = self._lower_kwargs
+        upper_kwargs_list = self._upper_kwargs
+        fixed_kwargs_list = self.fixed_kwargs_updated_list
+        for i in range(len(init_kwargs_list)):
+            init_kwargs, lower_kwargs, upper_kwargs, fixed_kwargs = (
+                init_kwargs_list[i],
+                lower_kwargs_list[i],
+                upper_kwargs_list[i],
+                fixed_kwargs_list[i],
+            )
+            if type(init_kwargs) is list:
+                # loop through dictionary
+                for k in range(len(init_kwargs)):
+                    _compare_bounds(
+                        init_kwargs[k],
+                        lower_kwargs[k],
+                        upper_kwargs[k],
+                        fixed_kwargs[k],
+                        model_index=k,
+                        model_class=_model_class_list[i],
+                    )
+            else:
+                _compare_bounds(
+                    init_kwargs,
+                    lower_kwargs,
+                    upper_kwargs,
+                    fixed_kwargs,
+                    model_index=None,
+                    model_class=_model_class_list[i],
+                )
+
     def set_init_state(self):
         """Set the current state of the parameters to the initial one.
 
@@ -673,3 +748,59 @@ class UpdateManager(object):
         :return: None
         """
         pass
+
+
+def _compare_bounds(
+    init_kwargs,
+    lower_kwargs,
+    upper_kwargs,
+    fixed_kwargs,
+    model_index=None,
+    model_class=None,
+):
+    """Raises Warnings when parameters are out of bounds.
+
+    :param init_kwargs: dictionary of initial parameters
+    :param lower_kwargs: lower bound dictionary
+    :param upper_kwargs: upper bound dictionary
+    :param model_index: integer of the index of the model in a certain model class
+    :param model_class: string of model class
+    :return: None
+    """
+    for key in init_kwargs:
+        # make sure linear parameters are not required when linear solver is on
+        if key not in fixed_kwargs:  # key not in ['amp'] or
+            if key not in lower_kwargs:
+                raise ValueError(
+                    "Variable %s missing in lower bounds of %s'th model of %s"
+                    % (key, model_index, model_class)
+                )
+            else:
+                if np.all(lower_kwargs[key] > init_kwargs[key]):
+                    Warning(
+                        "Variable %s of %s'th model of %s with initial guess %s is below the lower bound %s"
+                        % (
+                            key,
+                            model_index,
+                            model_class,
+                            init_kwargs[key],
+                            lower_kwargs[key],
+                        )
+                    )
+            if key not in upper_kwargs:
+                raise ValueError(
+                    "Variable %s missing in upper bounds of %s'th model of %s"
+                    % (key, model_index, model_class)
+                )
+            else:
+                if np.all(upper_kwargs[key] < init_kwargs[key]):
+                    Warning(
+                        "Variable %s of %s'th model of %s with initial guess %s is above the upper bound %s"
+                        % (
+                            key,
+                            model_index,
+                            model_class,
+                            init_kwargs[key],
+                            upper_kwargs[key],
+                        )
+                    )


Test patch:
diff --git a/test/test_Workflow/test_fitting_sequence.py b/test/test_Workflow/test_fitting_sequence.py
index 387efad7..932a956b 100644
--- a/test/test_Workflow/test_fitting_sequence.py
+++ b/test/test_Workflow/test_fitting_sequence.py
@@ -285,7 +285,7 @@ class TestFittingSequence(object):
             lens_temp[0]["theta_E"], self.kwargs_lens[0]["theta_E"], decimal=2
         )
 
-        logL = fittingSequence.best_fit_likelihood
+        logL = fittingSequence.best_fit_likelihood()
         print(logL, "test")
         # print(lens_temp, source_temp, lens_light_temp, ps_temp, special_temp)
         assert logL < 0
@@ -443,9 +443,9 @@ class TestFittingSequence(object):
         }
 
         lens_fixed = [{"center_x": 0.0, "center_y": 0.0}]
-        lens_sigma = [{"theta_E": 0.01}]
-        lens_lower = [{"theta_E": 0.1}]
-        lens_upper = [{"theta_E": 3.0}]
+        lens_sigma = [{"theta_E": 0.01, "center_x": 0.1, "center_y": 0.1}]
+        lens_lower = [{"theta_E": 0.1, "center_x": -10, "center_y": -10}]
+        lens_upper = [{"theta_E": 3.0, "center_x": 10, "center_y": 10}]
 
         source_fixed = [{}]
         source_sigma = [
diff --git a/test/test_Workflow/test_update_manager.py b/test/test_Workflow/test_update_manager.py
index 94d8ecb8..5e4ba7b9 100644
--- a/test/test_Workflow/test_update_manager.py
+++ b/test/test_Workflow/test_update_manager.py
@@ -209,6 +209,108 @@ class TestUpdateManager(object):
         self.manager.fix_image_parameters(image_index=0)
         assert 1 == 1
 
+    def test_check_initial_state(self):
+        kwargs_model = {
+            "lens_model_list": ["SHEAR"],
+        }
+        kwargs_constraints = {}
+        kwargs_likelihood = {}
+        kwargs_params = {}
+        lens_init = [{"e1": 0, "e2": 0}]
+        lens_sigma = [{"e1": 0.1, "e2": 0.1}]
+        lens_fixed = [{"ra_0": 0, "dec_0": 0}]
+        lens_lower = [{"e1": -1, "e2": -1}]
+        lens_upper = [{"e1": 1, "e2": 1}]
+        kwargs_params["lens_model"] = [
+            lens_init,
+            lens_sigma,
+            lens_fixed,
+            lens_lower,
+            lens_upper,
+        ]
+
+        manager = UpdateManager(
+            kwargs_model, kwargs_constraints, kwargs_likelihood, kwargs_params
+        )
+        # this should go smoothly through without warnings or raises
+        manager.check_initial_state()
+
+        # this should give a warning
+        lens_init = [{"e1": 0, "e2": 0}]
+        lens_sigma = [{"e1": 0.1, "e2": 0.1}]
+        lens_fixed = [{"ra_0": 0, "dec_0": 0}]
+        lens_lower = [
+            {"e1": 1, "e2": 1}
+        ]  # here we mix lower and upper bounds to cause a warning
+        lens_upper = [{"e1": -1, "e2": -1}]
+        kwargs_params["lens_model"] = [
+            lens_init,
+            lens_sigma,
+            lens_fixed,
+            lens_lower,
+            lens_upper,
+        ]
+
+        manager = UpdateManager(
+            kwargs_model, kwargs_constraints, kwargs_likelihood, kwargs_params
+        )
+        # should give a warning
+        manager.check_initial_state()
+
+        # should give a raise
+        # this should give a warning
+        lens_init = [{"e1": 0, "e2": 0}]
+        lens_sigma = [{"e1": 0.1, "e2": 0.1}]
+        lens_fixed = [{"ra_0": 0, "dec_0": 0}]
+        lens_lower = [
+            {
+                "e1": -1,
+            }
+        ]  # warning in lower bounds as keyword is missing
+        lens_upper = [{"e1": 1, "e2": 1}]
+        kwargs_params["lens_model"] = [
+            lens_init,
+            lens_sigma,
+            lens_fixed,
+            lens_lower,
+            lens_upper,
+        ]
+
+        manager = UpdateManager(
+            kwargs_model, kwargs_constraints, kwargs_likelihood, kwargs_params
+        )
+
+        from numpy.testing import assert_raises
+
+        with assert_raises(ValueError):
+            # should give a raise
+            manager.check_initial_state()
+
+        # should give a raise
+        # this should give a warning
+        lens_init = [{"e1": 0, "e2": 0}]
+        lens_sigma = [{"e1": 0.1, "e2": 0.1}]
+        lens_fixed = [{"ra_0": 0, "dec_0": 0}]
+        lens_lower = [{"e1": -1, "e2": -1}]
+        lens_upper = [{"e1": 1}]  # raise in upper bounds as keyword is missing
+        kwargs_params["lens_model"] = [
+            lens_init,
+            lens_sigma,
+            lens_fixed,
+            lens_lower,
+            lens_upper,
+        ]
+
+        manager = UpdateManager(
+            kwargs_model, kwargs_constraints, kwargs_likelihood, kwargs_params
+        )
+
+        from numpy.testing import assert_raises
+
+        with assert_raises(ValueError):
+            # should give a raise
+            manager.check_initial_state()
+
 
 if __name__ == "__main__":
     pytest.main()