Add bin/fg21sim and some updates to galactic/synchrotron

* Add new executable "bin/fg21sim"
* galactic/synchrotron: update to use "configs.get_path()"
* galactic/synchrotron: create output dir if not exists
* galactic/synchrotron: add logging support
* galactic/synchrotron: append FITS extension to filename
* galactic/synchrotron: pass the basic test

TODO:
* "output()" needs fixes with the FITS header

4 files changed, 102 insertions, 4 deletions

diff --git a/bin/fg21sim b/bin/fg21sim
new file mode 100755
index 0000000..8098ac3
--- /dev/null
+++ b/bin/fg21sim
@@ -0,0 +1,67 @@
+#!/usr/bin/env python3
+# -*- mode: python -*-
+#
+# Copyright (c) 2016 Weitian LI <liweitianux@live.com>
+# MIT license
+
+"""
+Simulate the low-frequency radio foregrounds for the 21cm EoR signal.
+"""
+
+import os
+import sys
+import argparse
+import logging
+
+import numpy as np
+
+from fg21sim.configs import configs, validate_configs
+from fg21sim.utils import setup_logging
+from fg21sim.galactic import Synchrotron as GalacticSynchrotron
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Simulate the radio foregrounds for 21cm EoR signal")
+    parser.add_argument("config", help="user configuration file")
+    parser.add_argument("-l", "--log", dest="loglevel", default=None,
+                        choices=["DEBUG", "INFO", "WARNING",
+                                 "ERROR", "CRITICAL"],
+                        help="set the log level")
+    parser.add_argument("-L", "--logfile", default=None,
+                        help="filename where to save the log messages")
+    parser.add_argument("-Q", "--quiet", action="store_true",
+                        help="be quiet so do not log messages to screen")
+    args = parser.parse_args()
+
+    configs.read_userconfig(args.config)
+    validate_configs(configs)
+
+    if args.quiet:
+        log_stream = ""
+    else:
+        log_stream = None
+
+    setup_logging(dict_config=configs.logging,
+                  level=args.loglevel,
+                  stream=log_stream,
+                  logfile=args.logfile)
+    tool = os.path.basename(sys.argv[0])
+    logger = logging.getLogger(tool)
+    logger.info("COMMAND: {0}".format(" ".join(sys.argv)))
+
+    frequencies = np.array(configs.frequencies, ndmin=1)
+    freq_unit = configs.getn("frequency/unit")
+    logger.info("Simulation frequencies: "
+                "{min:.2f} - {max:.2f} {unit} (#{num:d})".format(
+                    min=frequencies.min(), max=frequencies.max(),
+                    num=len(frequencies), unit=freq_unit))
+
+    logger.info("Simulating the Galactic synchrotron foreground ...")
+    gsynchrotron = GalacticSynchrotron(configs)
+    map_gsync = gsynchrotron.simulate(frequencies)
+    logger.info("Done simulate Galactic synchrotron foreground!")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/fg21sim/galactic/__init__.py b/fg21sim/galactic/__init__.py
new file mode 100644
index 0000000..4fe2659
--- /dev/null
+++ b/fg21sim/galactic/__init__.py
@@ -0,0 +1,4 @@
+# Copyright (c) 2016 Weitian LI <liweitianux@live.com>
+# MIT license
+
+from .synchrotron import Synchrotron
diff --git a/fg21sim/galactic/synchrotron.py b/fg21sim/galactic/synchrotron.py
index d05292f..614aae9 100644
--- a/fg21sim/galactic/synchrotron.py
+++ b/fg21sim/galactic/synchrotron.py
@@ -6,6 +6,7 @@ Diffuse Galactic synchrotron emission (unpolarized) simulations.
 """
 
 import os
+import logging
 from datetime import datetime, timezone
 
 import numpy as np
@@ -14,6 +15,9 @@ import astropy.units as au
 import healpy as hp
 
 
+logger = logging.getLogger(__name__)
+
+
 class Synchrotron:
     """Simulate the diffuse Galactic synchrotron emission based on an
     existing template.
@@ -40,20 +44,21 @@ class Synchrotron:
 
     def _set_configs(self):
         """Load the configs and set the corresponding class attributes."""
-        self.template_path = self.configs.getn(
+        self.template_path = self.configs.get_path(
             "galactic/synchrotron/template")
         self.template_freq = self.configs.getn(
             "galactic/synchrotron/template_freq")
         self.template_unit = au.Unit(
             self.configs.getn("galactic/synchrotron/template_unit"))
-        self.indexmap_path = self.configs.getn(
+        self.indexmap_path = self.configs.get_path(
             "galactic/synchrotron/indexmap")
         self.smallscales = self.configs.getn(
             "galactic/synchrotron/add_smallscales")
         # output
         self.prefix = self.configs.getn("galactic/synchrotron/prefix")
         self.save = self.configs.getn("galactic/synchrotron/save")
-        self.output_dir = self.configs.getn("galactic/synchrotron/output_dir")
+        self.output_dir = self.configs.get_path(
+            "galactic/synchrotron/output_dir")
         self.filename_pattern = self.configs.getn("output/filename_pattern")
         self.use_float = self.configs.getn("output/use_float")
         self.clobber = self.configs.getn("output/clobber")
@@ -63,6 +68,8 @@ class Synchrotron:
         self.lmax = self.configs.getn("common/lmax")
         # unit of the frequency
         self.freq_unit = au.Unit(self.configs.getn("frequency/unit"))
+        #
+        logger.info("Loaded and setup configurations")
 
     def _load_template(self):
         """Load the template map"""
@@ -71,6 +78,8 @@ class Synchrotron:
         self.template_header = fits.header.Header(header)
         self.template_nside = self.template_header["NSIDE"]
         self.template_ordering = self.template_header["ORDERING"]
+        logger.info("Loaded template map from {0} (Nside={1})".format(
+            self.template_path, self.template_nside))
 
     def _load_indexmap(self):
         """Load the spectral index map"""
@@ -79,6 +88,8 @@ class Synchrotron:
         self.indexmap_header = fits.header.Header(header)
         self.indexmap_nside = self.indexmap_header["NSIDE"]
         self.indexmap_ordering = self.indexmap_header["ORDERING"]
+        logger.info("Loaded spectral index map from {0} (Nside={1})".format(
+            self.indexmap_path, self.indexmap_nside))
 
     def _process_template(self):
         """Upgrade/downgrade the template to have the same Nside as
@@ -88,6 +99,8 @@ class Synchrotron:
         else:
             # upgrade/downgrade the resolution
             self.hpmap = hp.ud_grade(self.template, nside_out=self.nside)
+            logger.info("Upgrade/downgrade template map from Nside "
+                        "{0} to {1}".format(self.template_nside, self.nside))
 
     def _process_indexmap(self):
         """Upgrade/downgrade the spectral index map to have the same
@@ -97,6 +110,8 @@ class Synchrotron:
         else:
             # upgrade/downgrade the resolution
             self.hpmap_index = hp.ud_grade(self.indexmap, nside_out=self.nside)
+            logger.info("Upgrade/downgrade spectral index map from Nside "
+                        "{0} to {1}".format(self.indexmap_nside, self.nside))
 
     def _add_smallscales(self):
         """Add fluctuations on small scales to the template map.
@@ -128,11 +143,12 @@ class Synchrotron:
                        1.0 - np.exp(-ell[ell_idx]**2 * sigma_temp**2))
         cl[ell < self.lmin] = cl[self.lmin]
         # generate a realization of the Gaussian random field
-        gss = hp.synfast(cls=cl, nside=self.nside)
+        gss = hp.synfast(cls=cl, nside=self.nside, new=True)
         # whiten the Gaussian random field
         gss = (gss - gss.mean()) / gss.std()
         self.hpmap_smallscales = alpha * gss * self.hpmap**beta
         self.hpmap += self.hpmap_smallscales
+        logger.info("Added small-scale fluctuations")
 
     def _transform_frequency(self, frequency):
         """Transform the template map to the requested frequency,
@@ -149,6 +165,7 @@ class Synchrotron:
         header = fits.Header()
         header["COMP"] = ("Galactic synchrotron (unpolarized)",
                           "Emission component")
+        header["CREATOR"] = (__name__, "File creator")
         # TODO:
         history = []
         comments = []
@@ -157,6 +174,7 @@ class Synchrotron:
         for cmt in comments:
             header.add_comment(cmt)
         self.header = header
+        logger.info("Created FITS header")
 
     def output(self, hpmap, frequency):
         """Write the simulated synchrotron map to disk with proper
@@ -167,8 +185,13 @@ class Synchrotron:
             np.dtype("float64"): "D",
         }
         #
+        if not os.path.exists(self.output_dir):
+            os.mkdir(self.output_dir)
+            logger.info("Created output dir: {0}".format(self.output_dir))
+        #
         filename = self.filename_pattern.format(prefix=self.prefix,
                                                 frequency=frequency)
+        filename += ".fits"
         filepath = os.path.join(self.output_dir, filename)
         if not hasattr(self, "header"):
             self._make_header()
@@ -185,6 +208,7 @@ class Synchrotron:
                         format=FITS_COLUMN_FORMATS.get(hpmap.dtype))
         ], header=header)
         hdu.writeto(filepath, clobber=self.clobber, checksum=True)
+        logger.info("Write simulated map to file: {0}".format(filepath))
 
     def simulate(self, frequencies):
         """Simulate the synchrotron map at the specified frequencies."""
@@ -196,6 +220,8 @@ class Synchrotron:
         #
         hpmaps = []
         for f in np.array(frequencies, ndmin=1):
+            logger.info("Simulating synchrotron map at {0} ({1}) ...".format(
+                f, self.freq_unit))
             hpmap_f = self._transform_frequency(f)
             hpmaps.append(hpmap_f)
             if self.save:
diff --git a/setup.py b/setup.py
index e3d8c03..3d9a6d6 100755
--- a/setup.py
+++ b/setup.py
@@ -45,6 +45,7 @@ setup(
     ],
     packages=find_packages(exclude=["docs", "tests"]),
     scripts=[
+        "bin/fg21sim",
         "bin/healpix2hpx",
         "bin/hpx2healpix",
     ],