Support for unspecified properties in Observation and Instrument

coolest/api/composable_models.py

@@ -363,12 +363,26 @@ def evaluate_deflection(self, x, y):
             alpha_y += a_y
         return alpha_x, alpha_y
 
-    def evaluate_convergence(self, x, y):
+    def evaluate_convergence(self, x, y, mode='point', last_n_samples=None):
         """Evaluates the lensing convergence (i.e., 2D mass density) at given coordinates"""
-        kappa = np.zeros_like(x)
-        for k, (profile, params) in enumerate(zip(self.profile_list, self.param_list)):
-            kappa += profile.convergence(x, y, **params)
-        return kappa
+        self._check_eval_mode(mode)
+        if mode == 'point' or self._posterior_bool is False:
+            return self._eval_conv_point(x, y, self.param_list)
+        elif mode == 'posterior':
+            return self._eval_conv_posterior(x, y, self.post_param_list, last_n_samples)
+
+    def _eval_conv_point(self, x, y, param_list):
+        psi = np.zeros_like(x)
+        for k, profile in enumerate(self.profile_list):
+            psi += profile.convergence(x, y, **param_list[k])
+        return psi
+
+    def _eval_conv_posterior(self, x, y, param_list, last_n_samples):
+        # map the point function at each sample
+        use_all_samples = last_n_samples is None or last_n_samples <= 0
+        val_list = param_list if use_all_samples else param_list[-last_n_samples:]
+        mapped = map(partial(self._eval_conv_point, x, y), val_list)
+        return np.array(list(mapped))
 
     def evaluate_hessian(self, x, y):
         """Evaluates the lensing Hessian components at given coordinates"""
@@ -403,6 +417,20 @@ def ray_shooting(self, x, y):
         x_rs, y_rs = x - alpha_x, y - alpha_y
         return x_rs, y_rs
 
+    def shear(self, x, y):
+        """evaluates the reduce shear gamma"""
+        H_xx, H_xy, _, H_yy = self.evaluate_hessian(x, y)
+        gamma1 = 1.0 / 2 * (H_xx - H_yy)
+        gamma2 = H_xy
+        return gamma1, gamma2
+
+    def reduced_shear(self, x, y):
+        """evaluates the reduce shear g = gamma / 1 - kappa"""
+        kappa = self.evaluate_convergence(x, y)
+        gamma1, gamma2 = self.evaluate_convergence(x, y)
+        x_rs, y_rs = x - alpha_x, y - alpha_y
+        return gamma
+
 
 class ComposableLensModel(object):
     """Given a COOLEST object, evaluates a selection of entity and 

coolest/template/classes/grid.py

@@ -13,12 +13,10 @@
     'Grid',
     'PixelatedRegularGrid',
     'IrregularGrid',
+    'PixelatedRegularGridStack',
 ]
 
-SUPPORTED_CHOICES = [
-    'PixelatedRegularGrid',
-    'IrregularGrid',
-]
+SUPPORTED_CHOICES = list(set(__all__) - {'Grid'})
 
 
 class Grid(APIBaseObject):
@@ -295,3 +293,27 @@ def get_xyz(self, directory=None):
         y = data.field(1)
         z = data.field(2)
         return x, y, z
+
+
+class UnspecifiedGrid(Grid):
+    """Class that represents a generic grid of values, with just (optional)
+    FITS information. This may be used for e.g., as a workaround for storing 
+    other than standard photometric data, such as interferometric data.
+
+    Parameters
+    ----------
+    fits_path : str
+        Optional to some FITS file in which the values (and perhaps the coordinates)
+        are stored. This should be relative to the final COOLEST template file.
+    **kwargs_file : dic, optional
+        Any remaining keyword arguments for FitsFile
+    """
+
+    def __init__(self, 
+                 fits_path: str = None,
+                 **kwargs_file) -> None:
+        super().__init__(fits_path, **kwargs_file)
+        self.set_grid(
+            fits_path, 
+            check_fits_file=kwargs_file.get('check_fits_file', False)
+        )

coolest/template/classes/instrument.py

@@ -1,6 +1,6 @@
 __author__ = 'aymgal'
 
-from coolest.template.classes.psf import PSF
+from coolest.template.classes.psf import PSF, UnspecifiedPSF
 from coolest.template.classes.base import APIBaseObject
 
 
@@ -12,28 +12,29 @@ class Instrument(APIBaseObject):
     Parameters
     ----------
     pixel_size : float
-        Size in arcseconds of a single detector pixel.
+        Size in arcseconds of a single detector pixel, by default None.
     name : str, optional
         Name of the instrument, by default ""
     band : str, optional
         Name of the filter, by default ""
     readout_noise : float, optional
-        Readout noise (in electrons), by default 0.
+        Readout noise (in electrons) when it is relevant to the instrument, 
+        by default None.
     psf : PSF, optional
-        Instance of PSF object, by default None
+        Instance of PSF object, by default None (i.e. UnspecifiedPSF).
     """
 
     def __init__(self,
-                 pixel_size: float,
+                 pixel_size: float = None,
                  name: str = "",
                  band: str = "",
-                 readout_noise: float = 0.,
+                 readout_noise: float = None,
                  psf: PSF = None) -> None:
         self.name = name
         self.band = band
         self.pixel_size = pixel_size
         self.readout_noise = readout_noise
         if psf is None:
-            psf = PSF()
+            psf = UnspecifiedPSF()
         self.psf = psf
         super().__init__()

coolest/template/classes/noise.py

@@ -11,6 +11,7 @@
     'NoiseRealization',
     'InstrumentalNoise',
     'DrizzledNoise',
+    'UnspecifiedNoise',
 ]
 
 SUPPORTED_CHOICES = list(set(__all__) - {'Noise'})
@@ -124,3 +125,13 @@ def __init__(self, background_rms: float = 0.0,
         if wht_map is None:
             wht_map = PixelatedRegularGrid()
         super().__init__(ntype, background_rms=background_rms, wht_map=wht_map)
+
+
+class UnspecifiedNoise(Noise):
+    """Noise type that can be used when the noise can not be properly specified,
+    or noise properties are unknown. May be used for e.g. interferometric data.
+    """
+
+    def __init__(self) -> None:
+        ntype = self.__class__.__name__
+        super().__init__(ntype)

coolest/template/classes/observation.py

@@ -1,9 +1,10 @@
 __author__ = 'aymgal'
 
 from typing import Union
+import math
 
-from coolest.template.classes.grid import PixelatedRegularGrid
-from coolest.template.classes.noise import Noise
+from coolest.template.classes.grid import PixelatedRegularGrid, UnspecifiedGrid
+from coolest.template.classes.noise import Noise, UnspecifiedNoise
 from coolest.template.classes.base import APIBaseObject
 
 
@@ -37,20 +38,30 @@ def __init__(self,
                  # magnification_ratios: list = None
                  ) -> None:
         if pixels is None:
-            pixels = PixelatedRegularGrid()
+            pixels = UnspecifiedGrid()
         self.pixels = pixels
         self.exposure_time = exposure_time
         self.mag_zero_point = mag_zero_point          # magnitude zero-point (corresponds to 1 electron per second on the detector)
         self.mag_sky_brightness = mag_sky_brightness  # sky brightness (magnitude per arcsec^2)
         if noise is None:
-            noise = Noise()
+            noise = UnspecifiedNoise()
         self.noise = noise
         # self.time_delays = time_delays
         # self.magnification_ratios = magnification_ratios
         super().__init__()
 
     def check_consistency_with_instrument(self, instrument):
         """Checks that the data image is consistent with instrument properties"""
+        instru_pix_size = instrument.pixel_size
+        obs_pix_size = self.pixels.pixel_size
+        if instru_pix_size is None:
+            # when the instrument pixel size is undefined, we don't check anything
+            return
+        isclose_bool = math.isclose(instru_pix_size, obs_pix_size,
+                                    rel_tol=1e-09, abs_tol=0.0)
+        if obs_pix_size not in (0, None) and not isclose_bool:
+            raise ValueError(f"Pixel size of observation ({obs_pix_size}) is inconsistent with "
+                             f"the instrument pixel size ({instru_pix_size})")
         width  = abs(self.pixels.field_of_view_x[1] - self.pixels.field_of_view_x[0])
         height = abs(self.pixels.field_of_view_y[1] - self.pixels.field_of_view_y[0])
         num_pix_ra = round(width / instrument.pixel_size)

coolest/template/classes/psf.py

@@ -8,6 +8,7 @@
     'PSF',
     'PixelatedPSF',
     'GaussianPSF',
+    'UnspecifiedPSF',
 ]
 
 SUPPORTED_CHOICES = list(set(__all__) - {'PSF'})
@@ -75,3 +76,18 @@ def __init__(self, fwhm: float = 0.0,
                  description: str = None) -> None:
         psf_type = self.__class__.__name__
         super().__init__(psf_type, description=description, fwhm=fwhm)
+
+
+class UnspecifiedPSF(PSF):
+    """Can be used when the PSF is not, or can not be specified. May be used
+    with e.g. interferometric data.
+
+    Parameters
+    ----------
+    description : str, optional
+        Any details regarding the way the PSF has been estimated, by default None
+    """
+
+    def __init__(self, description: str = None) -> None:
+        psf_type = self.__class__.__name__
+        super().__init__(psf_type, description=description)

coolest/template/json.py

@@ -1,7 +1,6 @@
 import os
 import json
 import jsonpickle
-import math
 
 from coolest.template.standard import COOLEST
 from coolest.template.lazy import *
@@ -88,7 +87,7 @@ def dump_jsonpickle(self):
         with open(json_path, 'w') as f:
             f.write(result)
 
-    def load(self, skip_jsonpickle=False, validate=True, verbose=True):
+    def load(self, skip_jsonpickle=False, as_object=True, validate=True, verbose=True):
         """Read the JSON template file and build up the corresponding COOLEST object.
         It will first try to load the '_pyAPI' template if it exists using `jsonpickle`, 
         otherwise it will fall back to reading the pure json template.
@@ -97,6 +96,10 @@ def load(self, skip_jsonpickle=False, validate=True, verbose=True):
         ----------
         skip_jsonpickle : bool, optional
             If True, will not try to read the _pyAPI template with jsonpickle first, by default False
+        as_object : bool, optional
+            If False, returns the JSON as a plain dictionary instead, by default True
+        validate : bool, optional
+            If True, performs some consistency checks on the COOLEST object, by default True
         verbose : bool, optional
             If True, prints useful output for debugging, by default False
 
@@ -108,12 +111,15 @@ def load(self, skip_jsonpickle=False, validate=True, verbose=True):
         json_path = self.path + '.json'
         jsonpickle_path = self.path + self._api_suffix + '.json'
         if os.path.exists(jsonpickle_path) and not skip_jsonpickle:
+            if not as_object:
+                raise ValueError("Cannot load the JSON template with jsonpickle as a dictionary")
             instance = self.load_jsonpickle(jsonpickle_path)
         else:
             if verbose:
                 print(f"Template file '{jsonpickle_path}' not found, now trying to read '{json_path}'.")
-            instance = self.load_simple(json_path, as_object=True, validate=validate)
-        assert isinstance(instance, COOLEST)
+            instance = self.load_simple(json_path, as_object=as_object, validate=validate)
+        if as_object:
+            assert isinstance(instance, COOLEST)
         return instance
 
     def load_simple(self, json_path, as_object=True, validate=True):
@@ -124,7 +130,7 @@ def load_simple(self, json_path, as_object=True, validate=True):
         json_path: str
             Path to the json file to be read.
         as_object : bool, optional
-            _description_, by default True
+            If False, returns the JSON as a plain dictionary instead, by default True
 
         Returns
         -------
@@ -228,15 +234,6 @@ class constructors called during instantiation of the COOLEST object.
         ValueError
             In case observed instrumental pixel sizes are inconsistent
         """
-        # PIXEL SIZE
-        instru_pix_size = coolest.instrument.pixel_size
-        obs_pix_size = coolest.observation.pixels.pixel_size
-        isclose_bool = math.isclose(instru_pix_size, obs_pix_size,
-                                    rel_tol=1e-09, abs_tol=0.0)
-        if obs_pix_size not in (0, None) and not isclose_bool:
-            raise ValueError(f"Pixel size of observation ({obs_pix_size}) is inconsistent with "
-                             f"the instrument pixel size ({instru_pix_size})")
-        # INSTANCE METHODS
         coolest.observation.check_consistency_with_instrument(coolest.instrument)
         if coolest.likelihoods is not None:
             coolest.likelihoods.check_consistency_with_observation(coolest.observation)