Add a simple gaussian test

changes/165.phrosty.rst

@@ -0,0 +1 @@
+Add simple gaussian test

docs/development.rst

@@ -107,7 +107,7 @@ With these directories in place, run a container with::
 
 **on NERSC Perlmutter**, the command would be::
 
-  podman-hpc run -gpu -it \
+  podman-hpc run --gpu -it \
     --mount type=bind,source=$PWD,target=/home \
     --mount type=bind,source=$PWD/photometry_test_data,target=/photometry_test_data \
     --mount type=bind,source=$PWD/phrosty_temp,target=/phrosty_temp \
@@ -144,10 +144,12 @@ You should see all the options you can pass to phrosty.  There are a lot, becaus
 Try running::
 
   SNPIT_CONFIG=phrosty/tests/phrosty_test_config.yaml python phrosty/pipeline.py \
+    --oc ou2024 \
     --oid 20172782 \
     --ra 7.551093401915147 \
     --dec -44.80718106491529 \
     -b Y106 \
+    --ic ou2024 \
     -t phrosty/tests/20172782_instances_templates_1.csv \
     -s phrosty/tests/20172782_instances_science_2.csv \
     -p 3 -w 3 \

examples/perlmutter/phrosty_config.yaml

@@ -26,3 +26,4 @@ photometry:
 
     psf:
       type: ou24PSF
+      params: {}

phrosty/imagesubtraction.py

@@ -1,6 +1,7 @@
 __all__ = [ 'sky_subtract', 'stampmaker' ]
 
 # IMPORTS Standard:
+from astropy.io import fits
 import numpy as np
 import pathlib
 import random
@@ -63,11 +64,12 @@ def sky_subtract( img, temp_dir=None ):
         if isinstance( origimg, snappl.image.FITSImageOnDisk ):
             img = origimg.uncompressed_version( include=['data'] )
         else:
-            img = snappl.image.FITSImageOnDisk( path=tmpimpath )
+            img = snappl.image.FITSImage( path=tmpimpath, header=fits.header.Header() )
             img.data = origimg.data
-            img.save_data()
+            img.save( which='data' )
 
         SNLogger.debug( "Calling SEx_SkySubtract.SSS..." )
+        radius_cut_detmask = Config.get().value( 'photometry.phrosty.sfft.radius_cut_detmask' )
         ( _SKYDIP, _SKYPEAK, _PixA_skysub,
           _PixA_sky, PixA_skyrms ) = SEx_SkySubtract.SSS(FITS_obj=img.path,
                                                          FITS_skysub=tmpsubpath,
@@ -77,11 +79,12 @@ def sky_subtract( img, temp_dir=None ):
                                                          BACK_SIZE=64, BACK_FILTERSIZE=3,
                                                          DETECT_THRESH=1.5, DETECT_MINAREA=5,
                                                          DETECT_MAXAREA=0,
+                                                         RADIUS_CUT_DETMASK=radius_cut_detmask,
                                                          VERBOSE_LEVEL=2, MDIR=None)
         SNLogger.debug( "...back from SEx_SkySubtract.SSS" )
 
-        subim = snappl.image.FITSImageOnDisk( path=tmpsubpath )
-        detmaskim = snappl.image.FITSImageOnDisk( path=tmpdetmaskpath )
+        subim = snappl.image.FITSImage( path=tmpsubpath )
+        detmaskim = snappl.image.FITSImage( path=tmpdetmaskpath )
         skyrms = np.median( PixA_skyrms )
         return subim, detmaskim, skyrms
 
@@ -145,10 +148,15 @@ def stampmaker(ra, dec, shape, img, savedir=None, savename=None):
         if isinstance( origimg, snappl.image.FITSImageOnDisk ):
             img = origimg.uncompressed_version( include=['data'] )
         else:
-            barf = "".join( random.choices( "0123456789abcdef:", k=10 ) )
-            img = snappl.image.FITSImageOnDisk( path=savedir / f"{barf}.fits" )
+            barf = "".join( random.choices( "0123456789abcdef", k=10 ) )
+            # NOTE : this next line will break if not using an image type that
+            #   can return a FITS header!  The real solution is to fix SFFT
+            #   so that it's not dependent on FITS images; just pass what's
+            #   needed to Stamp_Generator.SG instead of assuming it will read
+            #   all the right things out of the header.
+            img = snappl.image.FITSImage( path=savedir / f"{barf}.fits", header=origimg.get_fits_header() )
             img.data = origimg.data
-            img.save_data()
+            img.save_data( which='data' )
 
         # TODO : if Stamp_Generator.SG can take a Path in FITS_StpLst, remove the str()
         Stamp_Generator.SG(FITS_obj=img.path, EXTINDEX=0, COORD=pxradec, COORD_TYPE='IMAGE',

phrosty/pipeline.py

@@ -126,10 +126,12 @@ def get_psf( self, ra, dec ):
 
         if self.psfobj is None:
             psftype = self.config.value( 'photometry.phrosty.psf.type' )
+            psfparams = self.config.value( 'photometry.phrosty.psf.params' )
             self.psfobj = PSF.get_psf_object( psftype, x=x, y=y,
                                               band=self.image.band,
                                               pointing=self.image.pointing,
-                                              sca=self.image.sca )
+                                              sca=self.image.sca,
+                                              **psfparams )
 
         stamp = self.psfobj.get_stamp( x, y )
         if self.keep_intermediate:
@@ -141,12 +143,13 @@ def get_psf( self, ra, dec ):
     def keep_psf_data( self, psf_data ):
         self.psf_data = psf_data
 
-    def free( self ):
+    def free( self, free_psf_data=False ):
         """Try to free memory.  More might be done here."""
         self.image.free()
         self.skysub_img.free()
         self.detmask_img.free()
-        self.psf_data = None
+        if free_psf_data:
+            self.psf_data = None
 
 
 class Pipeline:
@@ -240,7 +243,7 @@ def _read_csv( self, csvfile ):
             if not re.search( r"^\s*path\s+pointing\s+sca\s+mjd\s+band\s*$", hdrline ):
                 raise ValueError( f"First line of list file {csvfile} didn't match what was expected." )
             for line in ifp:
-                path, pointing, sca, mjd, band = line.split()
+                path, pointing, sca, _mjd, band = line.split()
                 if band == self.band:
                     # This should yell at us if the pointing or sca doesn't match what is read from the path
                     imlist.append( self.imgcol.get_image( path=path, pointing=pointing, sca=sca, band=band ) )
@@ -646,6 +649,11 @@ def __call__( self, through_step=None ):
         steps = steps[:stepdex+1]
 
         if 'sky_subtract' in steps:
+            # After this step is done, all images (both science and template)
+            #   will have the following fields set:
+            #     .skysub_img : sub-subtracted image
+            #     .detmask_img : deteciton mask image
+            #     .skyrms : float, median of sky image calculated by SExtractor
             SNLogger.info( "Running sky subtraction" )
             with nvtx.annotate( "skysub", color=0xff8888 ):
                 self.sky_sub_all_images()
@@ -654,6 +662,9 @@ def __call__( self, through_step=None ):
             SNLogger.info( f"After sky_subtract, memory usage = {tracemalloc.get_traced_memory()[1]/(1024**2):.2f} MB" )
 
         if 'get_psfs' in steps:
+            # After this step, all images (both science and template) will
+            #   have the .psf_data field set with the image-resolution PSF stamp data
+            #   (from a PSF.get_stamp() call.)
             SNLogger.info( "Getting PSFs" )
             with nvtx.annotate( "getpsfs", color=0xff8888 ):
                 self.get_psfs()
@@ -676,16 +687,63 @@ def log_fits_write_error( savepath, x ):
                     sfftifier = None
 
                     if 'align_and_preconvolve' in steps:
+                        # After this step, sfftifier will be a SpaceSFFT_CupyFlow
+                        #  object with the following fields filled.  All cupy arrays
+                        #  are float64 (even the DMASK!), and all of them are Transposed
+                        #  from what's stored in the PipelineImage object
+                        #    hdr_target : FITS header of science image
+                        #    hdr_object : FITS header of template image
+                        #    target_skyrms : float, median of sky for science image
+                        #    object_skyrms : float, median of sky for template image
+                        #    PixA_target_GPU  : science image data on GPU
+                        #    PixA_Ctarget_GPU : cross-convolved [with what?] science image on GPU
+                        #    PixA_object_GPU  : template image data on GPU
+                        #    PixA_resamp_object_GPU : warped template image data on GPU
+                        #    PixA_Cresampl_object_GPU : cross-convolved template image on GPU
+                        #    BlankMask_GPU : boolean array, True where PixA_resampl_object_GPU is 0.
+                        #    PixA_targetVar_GPU : noise CHECK THIS image for science image, on GPU
+                        #    PixA_objectVar_GPU : noise CHECK THIS image for template image, on GPU
+                        #    PiXA_resampl_objectVar_GPU : warped template variance data on GPU
+                        #    PixA_target_DMASK_GPU : detmask for science image, on GPU
+                        #    PixA_object_DMASK_GPU : detmask for template image, on GPU
+                        #    PixA_resampl_object_DMASK_GPU : warped dmask for template image, on GPU
+                        #    PSF_target_GPU : PSF stamp for science image
+                        #    PSF_object_GPU : PSF stamp for template image
+                        #    PSF_resampl_object_GPU : PSF stamp for template image rotated & resampled, on GPU
+                        #    sci_is_target : True
+                        #    GKerHW : 9 (int half-width of matching kernel, full width is 2*GKerHW+1)
+                        #    KerPolyOrder : config value of photometry.phrosty.kerpolyorder
+                        #    BGPolyOrder : 0
+                        #    ConstPhotRatio : True
+                        #    CUDA_DEVICE_4SUBTRACT : '0'
+                        #    GAIN : 1.0
+                        #    RANDOM_SEED : 10086
                         SNLogger.info( "...align_and_preconvolve" )
                         with nvtx.annotate( "align_and_pre_convolve", color=0x8888ff ):
                             sfftifier = self.align_and_pre_convolve( templ_image, sci_image )
 
                     if 'subtract' in steps:
+                        # After this step is done, two more fields in sfftifier are set:
+                        #    Solution_GPU  : --something--??
+                        #    PixA_DIFF_GPU : difference image ??
+                        # Get Lei to write docs on PureCupy_Customized_Packet.PCCP so
+                        #   we can figure out what these are
                         SNLogger.info( "...subtract" )
                         with nvtx.annotate( "subtraction", color=0x44ccff ):
                             sfftifier.sfft_subtraction()
 
                     if 'find_decorrelation' in steps:
+                        # This step does ...
+                        # After it's done, the following fields of sfftifier are set:
+                        #   Solution : CPU copy of Solution_GPU
+                        #   FKDECO_GPU : result of PureCupy_Decorrelation_Calculator.PCDC
+                        # Get Lei to write documentation on P:ureCupy_DeCorrelation_Calculator
+                        #   so we can figure out what this is, but I THINK it's a
+                        #   kernel that is used to convolve with things to "decorrelate".
+                        #   (From what?)
+                        # In addition the two local varaibles diff_var and diff_var_path are set.
+                        #   diff_var : variance in difference image (I THINK), on GPU
+                        #   diff_var_path : where we want to write diff_var in self.dia_out_dir
                         SNLogger.info( "...find_decorrelation" )
                         with nvtx.annotate( "find_decor", color=0xcc44ff ):
                             sfftifier.find_decorrelation()
@@ -697,6 +755,11 @@ def log_fits_write_error( savepath, x ):
                             diff_var_path = self.dia_out_dir / f"diff_var_{mess}"
 
                     if 'apply_decorrelation' in steps:
+                        # After this step is done, the following FITS files will be written in self.dia_out_dir:
+                        #   decorr_diff_* : difference image convolved with FKDECO_GPU
+                        #   diff_var_* : variance of difference image convolved with FKDECO_GPU
+                        #   decorr_zpt_path_* : preconvolved science image (PixA_Ctarget_GPU) convolved with FKDECO_GPU
+                        #   decor_psf_path_* : PSF stamp for science image, convolved with FKDECO_GPU
                         mess = f"{sci_image.image.name}-{templ_image.image.name}"
                         decorr_psf_path = self.dia_out_dir / f"decorr_psf_{mess}"
                         decorr_zptimg_path = self.dia_out_dir / f"decorr_zptimg_{mess}"

phrosty/tests/conftest.py

@@ -143,7 +143,7 @@ def one_science_image( scope="session" ):
     try:
         img = FITSImageOnDisk( path=('/photometry_test_data/ou2024/images/simple_model/'
                                      'Y106/35198/Roman_TDS_simple_model_Y106_35198_2.fits.gz'),
-                               imhdu=1,
+                               imagehdu=1,
                                pointing=35198,
                                sca=2 ).uncompressed_version()
         yield img

phrosty/tests/phrosty_test_config.yaml

@@ -21,14 +21,18 @@ photometry:
       A25ePSF_path: /snappl/snappl/tests/testdata/A25ePSF
       temp_dir: /phrosty_temp
 
+      simdex_server: https://roman-desc-simdex.rknop.net
+
   phrosty:
-    image_type: ou2024fits
     force_sky_subtract: true
     keep_intermediate: true
     remove_temp_dir: false
     mem_trace: false
     kerpolyorder: 2
 
+    sfft:
+      radius_cut_detmask: 2
+
     paths:
       scratch_dir: /phrosty_temp
       temp_dir: /phrosty_temp
@@ -37,3 +41,4 @@ photometry:
 
     psf:
       type: ou24PSF
+      params: {}

phrosty/tests/test_imagesubtraction.py

@@ -9,7 +9,6 @@
 from astropy.wcs import WCS
 
 from snappl.image import FITSImageOnDisk
-import phrosty
 import phrosty.imagesubtraction
 
 
@@ -27,7 +26,7 @@ def test_sky_subtract( dia_out_dir ):
         fits.writeto( in_path, imdata, header=hdr )
         img = FITSImageOnDisk( path=in_path )
 
-        subim, detmask, skymedrms = phrosty.imagesubtraction.sky_subtract( img, temp_dir=dia_out_dir )
+        subim, _detmask, skymedrms = phrosty.imagesubtraction.sky_subtract( img, temp_dir=dia_out_dir )
         assert skymedrms == pytest.approx( 10., abs=0.2 )
         assert subim.data.mean() == pytest.approx( 0., abs=5 * 10. / 512. )   # 3σ
         assert subim.data.std() == pytest.approx( 10., rel=0.05 )