Merge branch 'master' into 'current_stable_for_tests'

checksum fixing in header

See merge request !20

.gitignore

@@ -7,4 +7,6 @@ dist/*
 *interp.c
 !*libshao.so
 *.out
-pnodes
\ No newline at end of file
+pnodes
+.coverage
+coverage
\ No newline at end of file

ObservationSim/Instrument/FocalPlane.py

 import galsim
 import numpy as np
 
+
 class FocalPlane(object):
-    def __init__(self, config=None, chip_list=None, survey_type='Photometric', bad_chips=None):
+    def __init__(self, chip_list=None, survey_type='Photometric', bad_chips=None):
         """Get the focal plane layout
         """
         self.nchips = 42
@@ -14,7 +15,7 @@ class FocalPlane(object):
             self.bad_chips = bad_chips
             for chip_id in bad_chips:
                 self.ignore_chips.append(chip_id)
-        
+
         if chip_list is not None:
             for i in range(42):
                 if not (i+1 in chip_list):
@@ -39,25 +40,25 @@ class FocalPlane(object):
             for i in range(1, 31):
                 self.ignore_chips.append(i)
 
-        if config is not None:
-            self.nchip_x    = config["nchip_x"]
-            self.nchip_y    = config["nchip_y"]
-            self.npix_tot_x = config["npix_tot_x"]
-            self.npix_tot_y = config["npix_tot_y"]
-            self.npix_gap_x = config["npix_gap_x"]
-            self.npix_gap_y = config["npix_gap_y"]
-            if "chipLabelIDs" in config:
-                self.chipLabelIDs = config["chipLabelIDs"]
-            if "bad_chips" in config:
-                self.bad_chips = config["bad_chips"]
-        else:
-            self.nchip_x = 6
-            self.nchip_y = 5
-            self.npix_tot_x = 59516
-            self.npix_tot_y = 49752
-            self.npix_gap_x = (534, 1309)
-            self.npix_gap_y = 898
-        
+        # if config is not None:
+        #     self.nchip_x    = config["nchip_x"]
+        #     self.nchip_y    = config["nchip_y"]
+        #     self.npix_tot_x = config["npix_tot_x"]
+        #     self.npix_tot_y = config["npix_tot_y"]
+        #     self.npix_gap_x = config["npix_gap_x"]
+        #     self.npix_gap_y = config["npix_gap_y"]
+        #     if "chipLabelIDs" in config:
+        #         self.chipLabelIDs = config["chipLabelIDs"]
+        #     if "bad_chips" in config:
+        #         self.bad_chips = config["bad_chips"]
+        # else:
+        self.nchip_x = 6
+        self.nchip_y = 5
+        self.npix_tot_x = 59516
+        self.npix_tot_y = 49752
+        self.npix_gap_x = (534, 1309)
+        self.npix_gap_y = 898
+
         self._getCenter()
 
     def _getCenter(self):
@@ -65,7 +66,7 @@ class FocalPlane(object):
         self.cen_pix_y = 0
 
     def getChipLabel(self, chipID):
-        return str("0%d"%chipID)[-2:]
+        return str("0%d" % chipID)[-2:]
 
     def isBadChip(self, chipID):
         """Check if chip #(chipID) on the focal plane is bad or not
@@ -89,7 +90,8 @@ class FocalPlane(object):
             WCS of the focal plane
         """
         if logger is not None:
-            logger.info("    Construct the wcs of the entire image mosaic using Gnomonic/TAN projection")
+            logger.info(
+                "    Construct the wcs of the entire image mosaic using Gnomonic/TAN projection")
         if (xcen == None) or (ycen == None):
             xcen = self.cen_pix_x
             ycen = self.cen_pix_y
@@ -98,13 +100,14 @@ class FocalPlane(object):
         # dvdx =  -np.sin(img_rot.rad) * pix_scale
         # dvdy =  +np.cos(img_rot.rad) * pix_scale
 
-        dudx =  -np.cos(img_rot.rad) * pix_scale
-        dudy =  +np.sin(img_rot.rad) * pix_scale
-        dvdx =  -np.sin(img_rot.rad) * pix_scale
-        dvdy =  -np.cos(img_rot.rad) * pix_scale
-        
+        dudx = -np.cos(img_rot.rad) * pix_scale
+        dudy = +np.sin(img_rot.rad) * pix_scale
+        dvdx = -np.sin(img_rot.rad) * pix_scale
+        dvdy = -np.cos(img_rot.rad) * pix_scale
+
         moscen = galsim.PositionD(x=xcen, y=ycen)
-        sky_center = galsim.CelestialCoord(ra=ra*galsim.degrees, dec=dec*galsim.degrees)
+        sky_center = galsim.CelestialCoord(
+            ra=ra*galsim.degrees, dec=dec*galsim.degrees)
         affine = galsim.AffineTransform(dudx, dudy, dvdx, dvdy, origin=moscen)
         WCS = galsim.TanWCS(affine, sky_center, units=galsim.arcsec)
 
@@ -115,10 +118,10 @@ class FocalPlane(object):
         The sky coverage of an area
         """
         r2d = 180.0/np.pi
-        s1 = wcs.toWorld(galsim.PositionD(x0,y0))
-        s2 = wcs.toWorld(galsim.PositionD(x0,y1))
-        s3 = wcs.toWorld(galsim.PositionD(x1,y0))
-        s4 = wcs.toWorld(galsim.PositionD(x1,y1))
+        s1 = wcs.toWorld(galsim.PositionD(x0, y0))
+        s2 = wcs.toWorld(galsim.PositionD(x0, y1))
+        s3 = wcs.toWorld(galsim.PositionD(x1, y0))
+        s4 = wcs.toWorld(galsim.PositionD(x1, y1))
         ra = [s1.ra.rad*r2d, s2.ra.rad*r2d, s3.ra.rad*r2d, s4.ra.rad*r2d]
         dec = [s1.dec.rad*r2d, s2.dec.rad*r2d, s3.dec.rad*r2d, s4.dec.rad*r2d]
 

ObservationSim/sim_steps/add_objects.py

@@ -206,6 +206,12 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
         except Exception as e:
             traceback.print_exc()
             self.chip_output.Log_error(e)
+            self.chip_output.Log_error(
+                "pointing: #%d, chipID: %d" % (pointing.id, chip.chipID))
+            if obj.type == "galaxy":
+                self.chip_output.Log_error("obj id: %s" % (obj.param['id']))
+                self.chip_output.Log_error("    e1: %.5f\n    e2: %.5f\n    size: %f\n    bfrac: %f\n    detA: %f\n    g1: %.5f\n    g2: %.5f\n" % (
+                    obj.param['e1'], obj.param['e2'], obj.param['size'], obj.param['bfrac'], obj.param['detA'], obj.param['g1'], obj.param['g2']))
 
         # Unload SED:
         obj.unload_SED()

ObservationSim/sim_steps/readout_output.py

@@ -82,17 +82,16 @@ def quantization_and_output(self, chip, filt, tel, pointing, catalog, obs_param)
     
     
     hdu1 = fits.PrimaryHDU(header=self.h_prim)
-    hdu1.add_checksum()
-    hdu1.header.comments['CHECKSUM'] = 'HDU checksum'
-    hdu1.header.comments['DATASUM'] = 'data unit checksum'
-    self.updateHeaderInfo(header_flag='ext', keys = ['DATASECT'], values = [str(chip.img.array.shape[1])+'x'+str(chip.img.array.shape[0])])
+ 
+    self.updateHeaderInfo(header_flag='ext', keys = ['DATASECT'], values = [str(chip.img.array.shape[1]) + 'x' + str(chip.img.array.shape[0])])
     hdu2 = fits.ImageHDU(chip.img.array, header=self.h_ext)
-    hdu2.add_checksum()
-    hdu2.header.comments['XTENSION'] = 'extension type'
-    hdu2.header.comments['CHECKSUM'] = 'HDU checksum'
-    hdu2.header.comments['DATASUM'] = 'data unit checksum'
     hdu2.header.comments["XTENSION"] = "image extension"
-    hdu1 = fits.HDUList([hdu1, hdu2])
-    hdu1.writeto(fname, output_verify='ignore', overwrite=True)
+    
+    hdu = fits.HDUList([hdu1, hdu2])
+    hdu[0].add_datasum(when='data unit checksum')
+    hdu[0].add_checksum(when='HDU checksum', override_datasum=True)
+    hdu[1].add_datasum(when='data unit checksum')
+    hdu[1].add_checksum(when='HDU checksum', override_datasum=True)
+    hdu.writeto(fname, output_verify='ignore', overwrite=True)
 
     return chip, filt, tel, pointing

config/config_overall.yaml

@@ -11,7 +11,7 @@
 # can add some of the command-line arguments here as well;
 # ok to pass either way or both, as long as they are consistent
 work_dir: "/public/home/fangyuedong/project/workplace/"
-run_name: "magnification_test"
+run_name: "no_nonlinearity_test"
 
 # Project cycle and run counter are used to name the outputs
 project_cycle: 9
@@ -32,11 +32,13 @@ run_option:
 catalog_options:
   input_path:
     cat_dir: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/"
-    star_cat: "starcat/"
+    # star_cat: "starcat/"
+    star_cat: "starcat_C9/"
     galaxy_cat: "qsocat/cat2CSSTSim_bundle-50sqDeg/"
 
   SED_templates_path:
-    star_SED: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/SpecLib.hdf5"
+    # star_SED: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/SpecLib.hdf5"
+    star_SED: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/starcat_C9/"
     galaxy_SED: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/sedlibs/"
     AGN_SED: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/qsocat/qsosed/"
 
@@ -44,7 +46,7 @@ catalog_options:
   star_only: NO
 
   # Only simulate galaxies?
-  galaxy_only: YES
+  galaxy_only: NO
 
   # rotate galaxy ellipticity
   rotateEll: 0. # [degree]
@@ -57,7 +59,8 @@ obs_setting:
   # if you just want to run default pointing:
   # - pointing_dir: null
   # - pointing_file: null
-  pointing_file: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/pointing50_C9/pointing_50_1_n.dat"
+  # pointing_file: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/pointing50_C9/pointing_50_1_n.dat"
+  pointing_file: "/public/home/fangyuedong/project/unit_test_data/pointing_for_test.dat"
 
   obs_config_file: "/public/home/fangyuedong/project/csst_msc_sim/config/obs_config_SCI_WIDE_phot.yaml"
 
@@ -65,7 +68,7 @@ obs_setting:
   # - give a list of indexes of pointings: [ip_1, ip_2...]
   # - run all pointings: null
   # Note: only valid when a pointing list is specified
-  run_pointings: [0]
+  run_pointings: [1]
 
   # Whether to enable astrometric modeling
   enable_astrometric_model: True

config/obs_config_SCI_WIDE_phot.yaml

@@ -16,7 +16,7 @@ obs_id: "00000001" # this setting will only be used if pointing list file is not
 # Define list of chips
 # run_chips: [6,7,8,9,11,12,13,14,15,16,17,18,19,20,22,23,24,25]
 #run_chips: [1,2,3,4,5,10,21,26,27,28,29,30]
-run_chips: [8]
+run_chips: [9]
 
 # Define observation sequence
 call_sequence:

run_C6.pbs

@@ -2,9 +2,9 @@
 
 #SBATCH -J CSSTSim
 #SBATCH -N 1
-#SBATCH --ntasks-per-node=40
+#SBATCH --ntasks-per-node=6
 #SBATCH -p debug
-#SBATCH --mem=240G
+#SBATCH --mem=60G
 
 module load mpi/hpcx/2.4.1/gcc-7.3.1
 date
@@ -12,4 +12,4 @@ date
 #限定单节点任务数
 srun hostname -s | sort -n | awk -F"-" '{print $2}' | uniq > pnodes
 
-mpirun -mca pml ucx -x  UCX_NET_DEVICES=mlx5_0:1 -machinefile pnodes -np 40 --map-by node python3 /public/home/fangyuedong/project/csst_msc_sim/run_sim.py  --config_file config_overall.yaml --catalog C6_50sqdeg -c /public/home/fangyuedong/project/csst_msc_sim/config
\ No newline at end of file
+mpirun -mca pml ucx -x  UCX_NET_DEVICES=mlx5_0:1 -machinefile pnodes -np 6 --map-by node python3 /public/home/fangyuedong/project/csst_msc_sim/run_sim.py  --config_file config_overall.yaml --catalog C6_50sqdeg_ns -c /public/home/fangyuedong/project/csst_msc_sim/config
\ No newline at end of file

tests/.coverage

-}q
(U	collectorqUcoverage v3.6b3qUlinesq}u.
\ No newline at end of file

tests/test_focalplane.py

+import unittest
+import os
+import galsim
+from ObservationSim.Instrument import FocalPlane, Chip
+
+
+class TestFocalPlane(unittest.TestCase):
+    def __init__(self, methodName='runTest'):
+        super(TestFocalPlane, self).__init__(methodName)
+        self.dataPath = os.path.join(
+            os.getenv('UNIT_TEST_DATA_ROOT'), 'csst_msc_sim/csst_fz_msc')
+        self.focal_plane = FocalPlane(
+            chip_list=['8'])
+        self.assertTrue(self.focal_plane.cen_pix_x == 0)
+        self.assertTrue(self.focal_plane.cen_pix_y == 0)
+        test_focal_plane_phot = FocalPlane(survey_type='Photometric')
+        test_focal_plane_spec = FocalPlane(survey_type='Spectroscopic')
+        test_focal_plane_FGS = FocalPlane(survey_type='FGS')
+        test_focal_plane_bad_chips = FocalPlane(bad_chips=["1"])
+
+    def test_fp_method(self):
+        wcs = self.focal_plane.getTanWCS(
+            192.8595, 0., 0.*galsim.degrees, 0.0074)
+        sky_coverage = self.focal_plane.getSkyCoverage(
+            wcs, x0=-1, x1=0, y0=-1, y1=0)
+        print(sky_coverage.area())
+        self.assertTrue(abs(sky_coverage.area() - 0.0074**2/(3600.**2)) < 1e13)
+
+
+if __name__ == '__main_':
+    unittest.main()