Address miscellaneous TODO items

crispy/IFS.py

@@ -119,7 +119,8 @@ def polychromeIFS(par, inWavelist, inputcube,
         end=True)
 
     nframes = inputcube.data.shape[0]
-    allweights = None  # TODO, What is this variable? Add a helpful comment. 
+    # Variable for storing weight maps for each wavelength slice (currently unused)
+    allweights = None 
 
     if inputcube.data.shape[0] != len(wavelist):
         log.error('Number of wavelengths does not match the number of input cube slices')
@@ -171,17 +172,18 @@ def polychromeIFS(par, inWavelist, inputcube,
     if lam_arr is None:
         lam_arr = np.loadtxt(par.wavecalDir + "lamsol.dat")[:, 0]
 
-    hires_arrs = []  # TODO, Add a comment that describes what this variable is supposed to represent. Rename all occurence of this variable to 'high_res_arrays'  
+    # List to store high-resolution PSF arrays for each wavelength
+    high_res_arrays = []  
     if par.gaussian:
         for i in range(len(lam_arr)):
-            hiresarr = get_sim_hires(par, lam_arr[i])  # TODO rename all occurences of this variable to 'high_res_array'
-            hires_arrs += [hiresarr]
+            high_res_array = get_sim_hires(par, lam_arr[i])  # TODO rename all occurences of this variable to 'high_res_array'
+            high_res_arrays += [high_res_array]
         log.info('Creating Gaussian PSFLet templates')
         upsample = 10
     else:
         try:
             hires_list = np.sort(glob.glob(par.wavecalDir + 'hires_psflets_lam???.fits'))
-            hires_arrs = [pyf.getdata(filename) for filename in hires_list]
+            high_res_arrays = [pyf.getdata(filename) for filename in hires_list]
             log.info('Loaded PSFLet templates')
         except BaseException:
             log.error('Failed loading the PSFLet templates')
@@ -193,15 +195,16 @@ def polychromeIFS(par, inWavelist, inputcube,
     # Fill out the polyimage array
     if not parallel:
         for i in range(len(waveList)):
-            # TODO, add an informational comment about what imagePlaneRot represents. 
+            # Process and rotate the image plane for this wavelength slice, 
+            # then scale by wavelength bin width for flux conservation
             imagePlaneRot = (wavelist_endpts[i + 1] - wavelist_endpts[i]) * \
                 processImagePlane(par, interpolatedInputCube.data[i], noRot)
             inputCube += [imagePlaneRot]
             polyimage[i] = propagateLenslets(par,
                                              imagePlaneRot,
                                              wavelist_endpts[i],
                                              wavelist_endpts[i + 1],
-                                             hires_arrs,
+                                             high_res_arrays,
                                              lam_arr,
                                              upsample,
                                              nlam,
@@ -227,7 +230,7 @@ def polychromeIFS(par, inWavelist, inputcube,
                             imagePlaneRot,
                             wavelist_endpts[i],
                                wavelist_endpts[i + 1],
-                               hires_arrs,
+                               high_res_arrays,
                                lam_arr,
                                upsample,
                                nlam,

crispy/configs/HCIFS/params.py

@@ -95,7 +95,7 @@ def __init__(self, codeRoot='../'):
 
         self.BW = 0.18              # Spectral bandwidth
         self.npixperdlam = 2.0       # Number of pixels per spectral resolution element
-        self.nchanperspec_lstsq = 2.0  # Nspec per pixel for least squares
+        self.nchanperspec_lstsq = 2.0  # num_wavelengths per pixel for least squares
         # Spectral resolving power (extracted cubes have twice)
         self.R = 50
 
@@ -175,6 +175,6 @@ def makeHeader(self):
         self.hdr.append(
             ('RESLSTSQ',
              self.nchanperspec_lstsq,
-             'Nspec per Nyq. sample for lstsq extraction'),
+             'num_wavelengths per Nyq. sample for lstsq extraction'),
             end=True)
         self.hdr.append(('R', self.R, 'Spectral resolution'), end=True)

crispy/configs/PISCES/params.py

@@ -78,7 +78,7 @@ def __init__(self,codeRoot='../'):
         self.distortPISCES=False    # If True, use measured PISCES distortion/dispersion
         self.BW = 0.18              # Spectral bandwidth (if distortPISCES==False)
         self.npixperdlam = 2.0        # Number of pixels per spectral resolution element
-        self.nchanperspec_lstsq = 2.0 # Nspec per pixel for least squares 
+        self.nchanperspec_lstsq = 2.0 # num_wavelengths per pixel for least squares 
         self.R = 70                 # Spectral resolving power (extracted cubes have twice)
 
         # carry-over old parameter names

crispy/configs/WFIRST/params.py

@@ -90,7 +90,7 @@ def __init__(self, codeRoot='../'):
 
         self.BW = 0.18              # Spectral bandwidth
         self.npixperdlam = 2.0       # Number of pixels per spectral resolution element
-        self.nchanperspec_lstsq = 2.0  # Nspec per pixel for least squares
+        self.nchanperspec_lstsq = 2.0  # num_wavelengths per pixel for least squares
         # Spectral resolving power (extracted cubes have twice)
         self.R = 50
 
@@ -170,6 +170,6 @@ def makeHeader(self):
         self.hdr.append(
             ('RESLSTSQ',
              self.nchanperspec_lstsq,
-             'Nspec per Nyq. sample for lstsq extraction'),
+             'num_wavelengths per Nyq. sample for lstsq extraction'),
             end=True)
         self.hdr.append(('R', self.R, 'Spectral resolution'), end=True)

crispy/tools/lenslet.py

@@ -72,7 +72,7 @@ def propagateLenslets(
         imageplane,
         lam1,
         lam2,
-        hires_arrs=None,
+        high_res_arrays=None,
         lam_arr=None,
         upsample=3,
         nlam=10,
@@ -83,7 +83,7 @@ def propagateLenslets(
     Function propagateLenslets
 
     This is the main propagation function. It puts the PSFLets where they belong on the detector.
-    It uses template PSFLets given in hires_arrs, and can use also a pre-determined wavelength
+    It uses template PSFLets given in high_res_arrays, and can use also a pre-determined wavelength
     solution through the allcoef argument.
 
     Parameters
@@ -142,21 +142,21 @@ def propagateLenslets(
         # to do this later, saving a factor of a few in time.
         ################################################################
 
-        if (hires_arrs is None) or (lam_arr is None):
+        if (high_res_arrays is None) or (lam_arr is None):
             log.error('No template PSFLets given!')
             return
         else:
-            hires = np.zeros((hires_arrs[0].shape))
+            hires = np.zeros((high_res_arrays[0].shape))
             if lam <= np.amin(lam_arr):
-                hires[:] = hires_arrs[0]
+                hires[:] = high_res_arrays[0]
             elif lam >= np.amax(lam_arr):
-                hires[:] = hires_arrs[-1]
+                hires[:] = high_res_arrays[-1]
             else:
                 i1 = np.amax(np.arange(len(lam_arr))[np.where(lam > lam_arr)])
                 i2 = i1 + 1
-                hires = hires_arrs[i1] * \
+                hires = high_res_arrays[i1] * \
                     (lam - lam_arr[i1]) / (lam_arr[i2] - lam_arr[i1])
-                hires += hires_arrs[i2] * \
+                hires += high_res_arrays[i2] * \
                     (lam_arr[i2] - lam) / (lam_arr[i2] - lam_arr[i1])
 
             for i in range(hires.shape[0]):
@@ -244,23 +244,23 @@ def propagateLenslets(
                 x_hires = xcen[i] * 1. / image.shape[1]
                 y_hires = ycen[i] * 1. / image.shape[0]
 
-                x_hires = x_hires * hires_arrs[0].shape[1] - 0.5
-                y_hires = y_hires * hires_arrs[0].shape[0] - 0.5
+                x_hires = x_hires * high_res_arrays[0].shape[1] - 0.5
+                y_hires = y_hires * high_res_arrays[0].shape[0] - 0.5
 
                 totweight = 0
 
                 if x_hires <= 0:
                     i1 = i2 = 0
-                elif x_hires >= hires_arrs[0].shape[1] - 1:
-                    i1 = i2 = hires_arrs[0].shape[1] - 1
+                elif x_hires >= high_res_arrays[0].shape[1] - 1:
+                    i1 = i2 = high_res_arrays[0].shape[1] - 1
                 else:
                     i1 = int(x_hires)
                     i2 = i1 + 1
 
                 if y_hires < 0:
                     j1 = j2 = 0
-                elif y_hires >= hires_arrs[0].shape[0] - 1:
-                    j1 = j2 = hires_arrs[0].shape[0] - 1
+                elif y_hires >= high_res_arrays[0].shape[0] - 1:
+                    j1 = j2 = high_res_arrays[0].shape[0] - 1
                 else:
                     j1 = int(y_hires)
                     j2 = j1 + 1

crispy/tools/locate_psflets.py

@@ -129,29 +129,30 @@ def geninterparray(self, lam, allcoef, order=3):
             coef = np.linalg.lstsq(log_wavelength_powers, allcoef[:, i])[0]
             self.interp_arr[:, i] = coef
 
-    def return_locations_short(self, coef, xindx, yindx):
-        '''
-        Returns the x,y detector location of a given lenslet for a given polynomial fit
-
-        Parameters
-        ----------
-        coef: lists floats
-            Polynomial coefficients of fit for a single wavelength
-        xindx: int
-            X index of lenslet in lenslet array
-        yindx: int
-            Y index of lenslet in lenslet array
-
-        Returns
-        -------
-        interp_x: float
-            X coordinate on the detector
-        interp_y: float
-            Y coordinate on the detector
-        '''
-        # TODO, where does this 'coeforder' come from? Does this parent function actually get called from anywhere? Because clearly it's going to throw an error if it does.
-        interp_x, interp_y = transform(xindx, yindx, coeforder, coef) 
-        return interp_x, interp_y
+    # COMMENTED OUT: Function commented out to avoid confusion until we understand its purpose
+    # def return_locations_short(self, coef, xindx, yindx):
+    #     '''
+    #     Returns the x,y detector location of a given lenslet for a given polynomial fit
+
+    #     Parameters
+    #     ----------
+    #     coef: lists floats
+    #         Polynomial coefficients of fit for a single wavelength
+    #     xindx: int
+    #         X index of lenslet in lenslet array
+    #     yindx: int
+    #         Y index of lenslet in lenslet array
+
+    #     Returns
+    #     -------
+    #     interp_x: float
+    #         X coordinate on the detector
+    #     interp_y: float
+    #         Y coordinate on the detector
+    #     '''
+    #     # TODO, where does this 'coeforder' come from? Does this parent function actually get called from anywhere? Because clearly it's going to throw an error if it does.
+    #     interp_x, interp_y = transform(xindx, yindx, coeforder, coef) 
+    #     return interp_x, interp_y
 
     def return_res(self, lam, allcoef, xindx, yindx,
                    order=3, lam1=None, lam2=None):
@@ -782,28 +783,29 @@ def corrval(coef, x, y, input_image, order, trimfrac=0.1, show_plots=False):
     return score
 
 
-def corrvalsum(coef, x, y, filtered, order, trimfrac=0.1, gsize=2):
-    # TODO, is this function used anywhere in this repo? If not, comment it out. 
-    # TODO, add doscring
-    _x, _y = transform(x, y, order, coef)
-    ydim, xdim = filtered.shape
-    s = 0.0
-    ry = np.reshape(_y, -1)
-    rx = np.reshape(_x, -1)
-    for i in range(len(ry)):
-        yi = ry[i]
-        xi = rx[i]
-        xmin = int(xi) - gsize
-        xmax = xmin + 2 * gsize
-        ymin = int(yi) - gsize
-        ymax = ymin + 2 * gsize
-        if ymin > 2 * gsize and xmin > 2 * gsize and xmax < xdim - 2 * gsize and ymax < ydim - 2 * gsize:
-            # dx = xi - int(xi)
-            # dy = yi - int(yi)
-            # s+=np.sum(simplepsf(size=2*gsize,fwhm=fwhm,offx=dx,offy=dy)*filtered[ymin:ymax,xmin:xmax])
-            # s+=np.sum(gausspsf(size=2*gsize,fwhm=fwhm,offx=dx,offy=dy)*filtered[ymin:ymax,xmin:xmax])
-            s += np.sum(filtered[ymin:ymax, xmin:xmax])
-    return -s
+# COMMENTED OUT: Function unused in repository, commented out as requested
+# def corrvalsum(coef, x, y, filtered, order, trimfrac=0.1, gsize=2):
+#     # TODO, is this function used anywhere in this repo? If not, comment it out. 
+#     # TODO, add doscring
+#     _x, _y = transform(x, y, order, coef)
+#     ydim, xdim = filtered.shape
+#     s = 0.0
+#     ry = np.reshape(_y, -1)
+#     rx = np.reshape(_x, -1)
+#     for i in range(len(ry)):
+#         yi = ry[i]
+#         xi = rx[i]
+#         xmin = int(xi) - gsize
+#         xmax = xmin + 2 * gsize
+#         ymin = int(yi) - gsize
+#         ymax = ymin + 2 * gsize
+#         if ymin > 2 * gsize and xmin > 2 * gsize and xmax < xdim - 2 * gsize and ymax < ydim - 2 * gsize:
+#             # dx = xi - int(xi)
+#             # dy = yi - int(yi)
+#             # s+=np.sum(simplepsf(size=2*gsize,fwhm=fwhm,offx=dx,offy=dy)*filtered[ymin:ymax,xmin:xmax])
+#             # s+=np.sum(gausspsf(size=2*gsize,fwhm=fwhm,offx=dx,offy=dy)*filtered[ymin:ymax,xmin:xmax])
+#             s += np.sum(filtered[ymin:ymax, xmin:xmax])
+#     return -s
 
 
 def locatePSFlets(inImage, mask, polyorder=2, sig=0.7, coef=None, trimfrac=0.1,
@@ -831,7 +833,16 @@ def locatePSFlets(inImage, mask, polyorder=2, sig=0.7, coef=None, trimfrac=0.1,
         fraction of lenslet outliers (high & low
         combined) to trim in the minimization.  Default 0.1
         (5% trimmed on the high end, 5% on the low end)
-    # TODO, add argument definitions for the other arugments not listed here. 
+    mask: ndarray
+        Mask array for the image
+    phi: float
+        Rotation angle for the PSFlet grid. Default np.arctan2(1.926, -1)
+    scale: float
+        Scale factor for the PSFlet grid. Default 15.02
+    nlens: int
+        Number of lenslets. Default 108
+    finesearch: int
+        Fine search parameter. Default 3 
 
     Returns
     -------

crispy/tools/reduction.py

@@ -130,7 +130,7 @@ def testReduction(par, name, ifsimage):
     return cube
 
 
-def calculateWaveList(par, lam_list=None, Nspec=None, method='lstsq'):
+def calculateWaveList(par, lam_list=None, num_wavelengths=None, method='lstsq'):
     '''
     Computes the wavelength lists corresponding to the center and endpoints of each
     spectral bin. Wavelengths are separated by a constant value in log space. Number of
@@ -145,7 +145,7 @@ def calculateWaveList(par, lam_list=None, Nspec=None, method='lstsq'):
             calibration. Otherwise, we could decide to focus on a smaller/larger region of
             the spectrum to retrieve. The final processed cubes will have bins centered
             on lam_midpts
-    Nspec: int
+    num_wavelengths: int
             If specified, forces the number of bins in the final cube (uses np.linspace)
 
     Returns
@@ -160,30 +160,30 @@ def calculateWaveList(par, lam_list=None, Nspec=None, method='lstsq'):
         lamlist = np.loadtxt(par.wavecalDir + "lamsol.dat")[:, 0]
     else:
         lamlist = lam_list
-    if Nspec is None:
+    if num_wavelengths is None:
         if method == 'lstsq':
-            Nspec = int(
+            num_wavelengths = int(
                 np.log(
                     max(lamlist) /
                     min(lamlist)) *
                 par.R *
                 par.nchanperspec_lstsq +
                 1)
         else:
-            Nspec = int(
+            num_wavelengths = int(
                 np.log(
                     max(lamlist) /
                     min(lamlist)) *
                 par.R *
                 par.npixperdlam +
                 1)
-    log.info('Reduced cube will have %d wavelength bins' % (Nspec - 1))
-#     lam_endpts = np.linspace(min(lamlist), max(lamlist), Nspec)
+    log.info('Reduced cube will have %d wavelength bins' % (num_wavelengths - 1))
+#     lam_endpts = np.linspace(min(lamlist), max(lamlist), num_wavelengths)
 #     lam_midpts = (lam_endpts[1:]+lam_endpts[:-1])/2.
     loglam_endpts = np.linspace(
         np.log(
             min(lamlist)), np.log(
-            max(lamlist)), Nspec)
+            max(lamlist)), num_wavelengths)
     loglam_midpts = (loglam_endpts[1:] + loglam_endpts[:-1]) / 2
     lam_endpts = np.exp(loglam_endpts)
     lam_midpts = np.exp(loglam_midpts)
@@ -237,7 +237,7 @@ def lstsqExtract(par, name, ifsimage, smoothandmask=True, ivar=True, dy=3,
     yindx = polychromekey[2].data
     good = polychromekey[3].data
 
-    lam_midpts, lam_endpts = calculateWaveList(par, method='lstsq', Nspec=psflets.shape[0]+1)
+    lam_midpts, lam_endpts = calculateWaveList(par, method='lstsq', num_wavelengths=psflets.shape[0]+1)
 
     if fitbkgnd:
         n_add = 1
@@ -913,7 +913,7 @@ def intOptimalExtract(par, name, IFSimage, smoothandmask=True, sum=False):
     """
 
     loc = PSFLets(load=True, infiledir=par.wavecalDir)
-    #Nspec = int(par.BW*par.npixperdlam*par.R)
+    #num_wavelengths = int(par.BW*par.npixperdlam*par.R)
     lam_midpts, scratch = calculateWaveList(par, method='optext')
 
     datacube = fitspec_intpix_np(