Merge pull request #200 from mr-superonion/shear_obj_zshear

Add Shear Object enabling Redshift-dependent Shear

.github/workflows/shear_meas_tests.yml

@@ -30,21 +30,15 @@ jobs:
         shell: bash -l {0}
         run: |
           conda config --set always_yes yes
-          conda install -q mamba
 
-          mamba install -q stackvana=0
+          conda install -q mamba==1.5.1
 
+          mamba install -q --file requirements.txt
           mamba install -q \
             flake8 \
             pytest \
-            numpy \
-            galsim \
-            "numba!=0.54.0" \
-            ngmix \
-            lsstdesc-weaklensingdeblending \
             fitsio \
             meds \
-            hexalattice \
             ngmix
 
           pip install --no-deps -e .

.github/workflows/test.yml

@@ -29,20 +29,15 @@ jobs:
         shell: bash -l {0}
         run: |
           conda config --set always_yes yes
-          conda install -q mamba
 
-          mamba install -q stackvana=0
+          conda install -q mamba==1.5.1
 
+          mamba install -q --file requirements.txt
           mamba install -q \
             flake8 \
             pytest \
-            numpy \
-            galsim \
-            "numba!=0.54.0" \
-            ngmix \
-            lsstdesc-weaklensingdeblending \
             fitsio \
-            hexalattice
+            ngmix
 
           pip install --no-deps -e .
 

.gitignore

@@ -105,3 +105,4 @@ venv.bak/
 .mypy_cache/
 pytest_session.txt
 outputs
+.virtual_documents/

descwl_shear_sims/__init__.py

@@ -14,6 +14,7 @@
 from . import objlists
 from . import surveys
 from . import shifts
+from . import shear
 from . import constants
 from . import artifacts
 from . import masking

descwl_shear_sims/galaxies.py

@@ -203,7 +203,7 @@ def get_objlist(self, *, survey):
             objlist.append(self._get_galaxy(flux))
             shifts.append(galsim.PositionD(sarray['dx'][i], sarray['dy'][i]))
 
-        return objlist, shifts
+        return objlist, shifts, None
 
     def _get_galaxy(self, flux):
         """
@@ -633,13 +633,15 @@ def __init__(self, *, rng, coadd_dim, buff=0, layout='random'):
             nobj=nobj,
         )
 
+        # randomly sample from the catalog
         num = len(self)
         self.indices = self.rng.randint(
             0,
             self._wldeblend_cat.size,
             size=num,
         )
 
+        # do a random rotation for each galaxy
         self.angles = self.rng.uniform(low=0, high=360, size=num)
 
     def __len__(self):
@@ -672,11 +674,14 @@ def get_objlist(self, *, survey):
         sarray = self.shifts_array
         objlist = []
         shifts = []
+        redshifts = []
         for i in range(len(self)):
             objlist.append(self._get_galaxy(builder, band, i))
             shifts.append(galsim.PositionD(sarray['dx'][i], sarray['dy'][i]))
+            index = self.indices[i]
+            redshifts.append(self._wldeblend_cat[index]["redshift"])
 
-        return objlist, shifts
+        return objlist, shifts, redshifts
 
     def _get_galaxy(self, builder, band, i):
         """

descwl_shear_sims/objlists.py

@@ -20,7 +20,7 @@ def get_objlist(*, galaxy_catalog, survey, star_catalog=None, noise=None):
         objlist is a list of galsim GSObject with transformations applied. Shifts
         is an array with fields dx and dy for each object
     """
-    objlist, shifts = galaxy_catalog.get_objlist(survey=survey)
+    objlist, shifts, redshifts = galaxy_catalog.get_objlist(survey=survey)
 
     if star_catalog is not None:
         assert noise is not None
@@ -39,6 +39,7 @@ def get_objlist(*, galaxy_catalog, survey, star_catalog=None, noise=None):
     return {
         'objlist': objlist,
         'shifts': shifts,
+        'redshifts': redshifts,
         'star_objlist': sobjlist,
         'star_shifts': sshifts,
         'bright_objlist': bright_objlist,

descwl_shear_sims/shear.py

@@ -0,0 +1,89 @@
+import galsim
+import numpy as np
+
+
+def _ternary(n, n_bins):
+    """CREDIT: https://stackoverflow.com/questions/34559663/\
+        convert-decimal-to-ternarybase3-in-python"""
+    if n == 0:
+        return '0'
+    nums = []
+    while n:
+        n, r = divmod(n, 3)
+        nums.append(str(r))
+    return ''.join(reversed(nums)).zfill(n_bins)
+
+
+class ShearConstant(object):
+    """
+    Constant shear in the full exposure
+    Parameters
+    ----------
+    g1, g2:    Constant shear distortion
+    """
+    def __init__(self, g1, g2):
+        self.g1 = g1
+        self.g2 = g2
+        self.shear = galsim.Shear(g1=self.g1, g2=self.g2)
+        return
+
+    def get_shear(self, redshift=None, shift=None):
+        """
+        Returns
+        ---------
+        shear (galsim.Shear)        shear distortion on the galaxy
+        """
+        return self.shear
+
+
+class ShearRedshift(object):
+    """
+    Constant shear in each redshift slice
+    """
+    def __init__(self, z_bounds, mode, g_dist="g1", shear_value=0.02):
+        assert isinstance(mode, int), "mode must be an integer"
+        self.nz_bins = int(len(z_bounds) - 1)
+        # nz_bins is the number of redshift bins
+        # note that there are three options in each redshift bin
+        # 0: g=-0.02; 1: g=0.02; 2: g=0.00
+        # for example, number of redshift bins is 4, (nz_bins = [0., 0.5, 1.0,
+        # 1.5, 2.0]) if mode = 7 which in ternary is "0021" --- meaning that
+        # the shear is (-0.02, -0.02, 0.00, 0.02) in each bin, respectively.
+        self.code = _ternary(int(mode), self.nz_bins)
+        assert 0 <= int(mode) < 3 ** self.nz_bins, "mode code is too large"
+        # maybe we need it to be more flexible in the future
+        # but now we keep the linear spacing
+        self.z_bounds = z_bounds
+        self.g_dist = g_dist
+        self.shear_value = shear_value
+        self.shear_list = self.determine_shear_list(self.code)
+        return
+
+    def determine_shear_list(self, code):
+        values = [-self.shear_value, self.shear_value, 0.0]
+        shear_list = [values[int(i)] for i in code]
+        return shear_list
+
+    def _get_zshear(self, redshift):
+        bin_num = np.searchsorted(self.z_bounds, redshift, side="left") - 1
+        nz = len(self.z_bounds) - 1
+        if bin_num < nz and bin_num >= 0:
+            # if the redshift is within the boundaries of lower and uper limits
+            # we add shear
+            shear = self.shear_list[bin_num]
+        else:
+            # if not, we set shear to 0 and leave the galaxy image undistorted
+            shear = 0.0
+        return shear
+
+    def get_shear(self, redshift, shift=None):
+        shear = self._get_zshear(redshift)
+        if self.g_dist == 'g1':
+            gamma1, gamma2 = (shear, 0.)
+        elif self.g_dist == 'g2':
+            gamma1, gamma2 = (0., shear)
+        else:
+            raise ValueError("g_dist must be either 'g1' or 'g2'")
+
+        shear_obj = galsim.Shear(g1=gamma1, g2=gamma2)
+        return shear_obj

descwl_shear_sims/shifts.py

@@ -57,6 +57,7 @@ def get_shifts(
         elif layout == 'random':
             # area covered by objects
             if nobj is None:
+                # in units of square arcmin
                 area = ((coadd_dim - 2*buff)*SCALE/60)**2
                 nobj_mean = max(area * RANDOM_DENSITY, 1)
                 nobj = rng.poisson(nobj_mean)