improve the simulation code

bin/fpfs_sim.py

@@ -15,10 +15,12 @@
 #
 import os
 import gc
+import glob
 import fpfs
 import json
 import galsim
 import schwimmbad
+import numpy as np
 from argparse import ArgumentParser
 from configparser import ConfigParser
 
@@ -27,138 +29,130 @@ class Worker(object):
     def __init__(self, config_name):
         cparser = ConfigParser()
         cparser.read(config_name)
-        self.simname = cparser.get("simulation", "sim_name")
+        self.sim_method = cparser.get("simulation", "sim_method")
+        self.gal_type = cparser.get("simulation", "gal_type").lower()
         self.imgdir = cparser.get("simulation", "img_dir")
-        self.infname = cparser.get("simulation", "input_name")
+        self.nrot = cparser.getint("simulation", "nrot")
+        self.band_name = cparser.get("simulation", "band")
         self.scale = cparser.getfloat("survey", "pixel_scale")
         self.image_nx = cparser.getint("survey", "image_nx")
         self.image_ny = cparser.getint("survey", "image_ny")
         assert self.image_ny == self.image_nx, "'image_nx' must equals 'image_ny'!"
-        if "basic" in self.simname or "small" in self.simname:
-            assert self.image_nx % 256 == 0, "'image_nx' must be divisible by 256 ."
-        self.psfInt = None
-        self.outdir = os.path.join(self.imgdir, self.simname)
+        self.psf_obj = None
+        self.outdir = os.path.join(self.imgdir, self.sim_method)
         if not os.path.exists(self.outdir):
             os.makedirs(self.outdir, exist_ok=True)
-        if "galaxy" in self.simname:
-            assert (
-                "basic" in self.simname
-                or "small" in self.simname
-                or "cosmo" in self.simname
-            )
-            if cparser.has_option("survey", "psf_fwhm"):
-                seeing = cparser.getfloat("survey", "psf_fwhm")
-                self.prepare_psf(seeing, psf_type="moffat")
-                print("Using modelled Moffat PSF with seeing %.2f arcsec. " % seeing)
-            else:
-                if not cparser.has_option("survey", "psf_filename"):
-                    raise ValueError("Do not have survey-psf_file option")
-                else:
-                    self.psffname = cparser.get("survey", "psf_filename")
-                    print("Using PSF from input file. ")
-            glist = []
-            if cparser.getboolean("distortion", "test_g1"):
-                glist.append("g1")
-            if cparser.getboolean("distortion", "test_g2"):
-                glist.append("g2")
-            if len(glist) > 0:
-                zlist = json.loads(cparser.get("distortion", "shear_z_list"))
-                self.pendList = ["%s-%s" % (i1, i2) for i1 in glist for i2 in zlist]
+        assert (
+            self.sim_method in ["fft", "mc"]
+        )
+        assert (
+            self.gal_type in ["mixed", "sersic", "bulgedisk"]
+        )
+        if cparser.has_option("survey", "psf_fwhm"):
+            seeing = cparser.getfloat("survey", "psf_fwhm")
+            self.prepare_psf(seeing, psf_type="moffat")
+            print("Using modelled Moffat PSF with seeing %.2f arcsec. " % seeing)
+        else:
+            if not cparser.has_option("survey", "psf_filename"):
+                raise ValueError("Do not have survey-psf_file option")
             else:
-                # this is for non-distorted image simulation
-                self.pendList = ["g1-1111"]
-            print(
-                "We will test the following constant shear distortion setups %s. "
-                % self.pendList
-            )
-            self.add_halo = cparser.getboolean("distortion", "add_halo")
-            self.shear_value = cparser.getfloat("distortion", "shear_value")
-            if self.add_halo:
-                assert self.pendList == [
-                    "g1-1111"
-                ], "Do not support adding both \
-                        shear from halo and constant shear."
-        elif "noise" in self.simname:
-            self.pendList = [0]
+                self.psffname = cparser.get("survey", "psf_filename")
+                print("Using PSF from input file. ")
+        glist = []
+        if cparser.getboolean("distortion", "test_g1"):
+            glist.append("g1")
+        if cparser.getboolean("distortion", "test_g2"):
+            glist.append("g2")
+        if len(glist) > 0:
+            zlist = json.loads(cparser.get("distortion", "shear_z_list"))
+            self.pendList = ["%s-%s" % (i1, i2) for i1 in glist for i2 in zlist]
         else:
-            raise ValueError(
-                "Cannot setup the task for sim_name=%s!!\\ \
-                    Must contain 'galaxy' or 'noise'"
-                % self.simname
-            )
+            # this is for non-distorted image simulation
+            self.pendList = ["g1-2"]
+        print(
+            "We will test the following constant shear distortion setups %s. "
+            % self.pendList
+        )
+        self.shear_value = cparser.getfloat("distortion", "shear_value")
+        self.rot_list = [np.pi / self.nrot * i for i in range(self.nrot)]
         return
 
     def prepare_psf(self, seeing, psf_type):
         psffname = os.path.join(self.outdir, "psf-%d.fits" % (seeing * 100))
         if psf_type.lower() == "moffat":
-            self.psfInt = galsim.Moffat(
+            self.psf_obj = galsim.Moffat(
                 beta=3.5, fwhm=seeing, trunc=seeing * 4.0
             ).shear(e1=0.02, e2=-0.02)
         else:
             raise ValueError("Only support moffat PSF.")
-        psf_image = self.psfInt.drawImage(nx=45, ny=45, scale=self.scale)
+        psf_image = self.psf_obj.shift(
+            0.5 * self.scale, 0.5 * self.scale
+        ).drawImage(nx=64, ny=64, scale=self.scale)
         psf_image.write(psffname)
         return
 
-    def run(self, Id):
-        if "noise" in self.simname:
-            # do pure noise image simulation
-            fpfs.simutil.make_noise_sim(self.outdir, self.infname, Id, scale=self.scale)
-        else:
-            if self.psfInt is None:
-                if "%" in self.psffname:
-                    psffname = self.psffname % Id
-                else:
-                    psffname = self.psffname
-                assert os.path.isfile(psffname), "Cannot find input PSF file"
-                psf_image = galsim.fits.read(psffname)
-                self.psfInt = galsim.InterpolatedImage(
-                    psf_image, scale=self.scale, flux=1.0
+    def run(self, ifield):
+        print("start ID: %d" % (ifield))
+        if self.psf_obj is None:
+            if "%" in self.psffname:
+                psffname = self.psffname % ifield
+            else:
+                psffname = self.psffname
+            assert os.path.isfile(psffname), "Cannot find input PSF file"
+            psf_image = galsim.fits.read(psffname)
+            self.psf_obj = galsim.InterpolatedImage(
+                psf_image, scale=self.scale, flux=1.0
+            )
+            del psf_image
+        for pp in self.pendList:
+            # do basic stamp-like image simulation
+            nfiles = len(glob.glob("%s/image-%05d_%s_rot*_%s.fits" % (
+                self.outdir,
+                ifield,
+                pp,
+                self.band_name,
+            )))
+            if nfiles == self.nrot:
+                print("We already have all the output files for %s" % pp)
+                continue
+            sim_img = fpfs.simutil.make_isolate_sim(
+                sim_method="fft",
+                psf_obj=self.psf_obj,
+                gname=pp,
+                seed=ifield,
+                ny=self.image_ny,
+                nx=self.image_nx,
+                scale=self.scale,
+                do_shift=False,
+                shear_value=self.shear_value,
+                nrot=1,
+                rot2=self.rot_list,
+                gal_type=self.gal_type,
+            )
+            for irot in range(self.nrot):
+                gal_fname = "%s/image-%05d_%s_rot%d_%s.fits" % (
+                    self.outdir,
+                    ifield,
+                    pp,
+                    irot,
+                    self.band_name,
                 )
-                del psf_image
-            for pp in self.pendList:
-                if "basic" in self.simname or "small" in self.simname:
-                    # do basic stamp-like image simulation
-                    fpfs.simutil.make_basic_sim(
-                        self.outdir,
-                        self.psfInt,
-                        pp,
-                        Id,
-                        catname=self.infname,
-                        ny=self.image_ny,
-                        nx=self.image_nx,
-                        scale=self.scale,
-                        shear_value=self.shear_value,
-                    )
-                elif "cosmo" in self.simname:
-                    # do blended cosmo-like image simulation
-                    fpfs.simutil.make_cosmo_sim(
-                        self.outdir,
-                        self.psfInt,
-                        pp,
-                        Id,
-                        catname=self.infname,
-                        ny=self.image_ny,
-                        nx=self.image_nx,
-                        scale=self.scale,
-                        shear_value=self.shear_value,
-                    )
-        gc.collect()
-        print("finish ID: %d" % (Id))
+                fpfs.io.save_image(gal_fname, sim_img[irot])
+            gc.collect()
+        print("finish ID: %d" % (ifield))
         return
 
-    def __call__(self, Id):
-        print("start ID: %d" % (Id))
-        return self.run(Id)
+    def __call__(self, ifield):
+        return self.run(ifield)
 
 
 if __name__ == "__main__":
     parser = ArgumentParser(description="fpfs simulation")
     parser.add_argument(
-        "--minId", required=True, type=int, help="minimum id number, e.g. 0"
+        "--min_id", required=True, type=int, help="minimum id number, e.g. 0"
     )
     parser.add_argument(
-        "--maxId", required=True, type=int, help="maximum id number, e.g. 4000"
+        "--max_id", required=True, type=int, help="maximum id number, e.g. 4000"
     )
     parser.add_argument("--config", required=True, type=str, help="configure file name")
     group = parser.add_mutually_exclusive_group()
@@ -177,7 +171,7 @@ def __call__(self, Id):
     pool = schwimmbad.choose_pool(mpi=args.mpi, processes=args.n_cores)
 
     worker = Worker(args.config)
-    refs = list(range(args.minId, args.maxId))
+    refs = list(range(args.min_id, args.max_id))
     # worker(1)
     for r in pool.map(worker, refs):
         pass

bin/configure_files/desc_sim_noise1.ini → bin/tests/config_n0.ini

@@ -1,8 +1,9 @@
 [procsim]
-img_dir     =   /lustre/work/xiangchong.li/work/FPFS2/sim_desc/basic
-psf_fname   =   /lustre/work/xiangchong.li/work/FPFS2/sim_desc/PSF_basic_32.fits
-cat_dir     =   cat_n1_p1
-sum_dir     =   sum_n1_p1
+img_dir     =   ./fft
+psf_fname   =   ./fft/psf-60.fits
+cat_dir     =   cat_n0
+sum_dir     =   sum_n0
+
 
 
 [FPFS]
@@ -15,6 +16,8 @@ alpha       =   0.35
 beta        =   0.92
 sigma_as    =   0.52
 sigma_det   =   0.53
+ncov_fname  =   cat_n1/cov_matrix.fits
+noise_rev   =   False
 
 
 [distortion]
@@ -24,7 +27,8 @@ shear_value =   0.02
 
 
 [survey]
-band = a
-mag_zero   =   30.
-noise_ratio =  1
-pixel_scale =  0.2
+band = i
+mag_zero   =   27.
+noise_std =    0.
+pixel_scale =  0.168
+

bin/tests/config_n1.ini

@@ -0,0 +1,33 @@
+[procsim]
+img_dir     =   ./fft
+psf_fname   =   ./fft/psf-60.fits
+cat_dir     =   cat_n1
+sum_dir     =   sum_n1
+
+
+
+[FPFS]
+nnord       =   4
+rcut        =   32
+ratio       =   1.596
+c0          =   2.46
+c2          =   22.74
+alpha       =   0.35
+beta        =   0.92
+sigma_as    =   0.52
+sigma_det   =   0.53
+noise_rev   =   True
+
+
+[distortion]
+test_g1     =   True
+test_g2     =   False
+shear_value =   0.02
+
+
+[survey]
+band = i
+mag_zero   =   27.
+noise_std =    3e-3
+pixel_scale =  0.168
+

bin/tests/config_sim_gal.ini

@@ -0,0 +1,27 @@
+[simulation]
+; type of simulation
+; stamp-based simulation
+sim_method  =   fft
+gal_type    =   Mixed
+img_dir     =   ./
+nrot        =   2
+band        =   i
+
+
+[distortion]
+add_halo    =   False
+test_g1     =   True
+test_g2     =   False
+; 0000 for negative distortion; 1111 for no distortion; and 2222 for positive
+; distortion
+shear_z_list=   ["0", "1"]
+shear_value =   0.02
+
+
+[survey]
+; survey info for HSC coadds
+pixel_scale =   0.168
+psf_fwhm    =   0.60
+mag_zero    =   27
+image_nx    =   1024
+image_ny    =   1024

docs/examples/notebooks/1_simulate_noise_image.ipynb → docs/examples/1_simulate_noise_image.ipynb

@@ -17,7 +17,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 2,
    "id": "51c95e6d-238c-4dd8-ba70-8123fbf44e9b",
    "metadata": {},
    "outputs": [
@@ -40,7 +40,7 @@
        "         0.0438892 , -0.1618046 ]], dtype=float32), wcs=galsim.OffsetWCS(0.168, galsim.PositionI(x=-1, y=-1), galsim.PositionD(x=0.0, y=0.0)))"
       ]
      },
-     "execution_count": 4,
+     "execution_count": 2,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -84,7 +84,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.6"
+   "version": "3.11.5"
   }
  },
  "nbformat": 4,