Merge branch 'ephem-handling' into 'master'

ephemerides handling fixes

See merge request GregAshton/PyFstat!25

README.md

@@ -12,9 +12,10 @@ we have a number of scripts demonstrating different use cases.
 
 ## Installation
 
-### `python` installation
-The scripts are written in `python 2.7+` and therefore require a working
-`python` installation. While many systems come with a system wide python
+### python installation
+This package works best with `python3.5+`,
+with higher versions to be required soon.
+While many systems come with a system wide python
 installation, it can often be easier to manage a user-specific python
 installation. This way one does not require root access to install or remove
 modules. One method to do this, is to use the `conda` system, either through
@@ -82,7 +83,7 @@ $ ./configure --prefix=${HOME}/lalsuite-install --disable-all-lal --enable-lalpu
 ```
 
 
-### `pyfstat` installation
+### PyFstat installation
 
 The module and associated scripts can be installed system wide (or to the currently active venv),
 assuming you are in the source directory, via
@@ -97,14 +98,37 @@ $ python -c 'import pyfstat'
 if no error message is output, then you have installed `pyfstat`. Note that
 the module will be installed to whichever python executable you call it from.
 
-### Ephemeris installation
-
-The scripts require paths to earth and sun ephemeris files in order to use the
-`lalpulsar.ComputeFstat` module. This should be automatically picked up from
-the $LALPULSAR_DATADIR environment variable, defaulting to the
-00-40-DE421 ephemerides or 00-19-DE421 as a backup.
-Alternatively, these can either be manually specified when initialising
-each search (as one of the arguments), or simply by placing a file
+### Ephemerides installation
+
+PyFstat requires paths to earth and sun ephemerides files
+in order to use the `lalpulsar.ComputeFstat` module and various `lalapps` tools.
+
+If you have done `pip install lalsuite`,
+you need to manually download at least these two files:
+*  [earth00-40-DE405.dat.gz](https://git.ligo.org/lscsoft/lalsuite/raw/master/lalpulsar/src/earth00-40-DE405.dat.gz)
+*  [sun00-40-DE405.dat.gz](https://git.ligo.org/lscsoft/lalsuite/raw/master/lalpulsar/src/sun00-40-DE405.dat.gz)
+
+(Other ephemerides versions exist, but these should be sufficient for most applications.)
+You then need to tell PyFstat where to find these files,
+by either setting an environment variable $LALPULSAR_DATADIR
+or by creating a `~/.pyfstat.conf` file as described further below.
+If you are working with a virtual environment,
+you should be able to get a full working ephemerides installation with these commands:
+```
+mkdir $VIRTUAL_ENV/share/lalpulsar
+wget https://git.ligo.org/lscsoft/lalsuite/raw/master/lalpulsar/src/earth00-40-DE405.dat.gz -P $VIRTUAL_ENV/share/lalpulsar
+wget https://git.ligo.org/lscsoft/lalsuite/raw/master/lalpulsar/src/sun00-40-DE405.dat.gz -P $VIRTUAL_ENV/share/lalpulsar
+echo 'export LALPULSAR_DATADIR=$VIRTUAL_ENV/share/lalpulsar' >> ${VIRTUAL_ENV}/bin/activate
+deactivate
+source path/to/venv/bin/activate
+```
+
+If instead you have built and installed lalsuite from source,
+and set your path up properly through something like
+`source $MYLALPATH/etc/lalsuite-user-env.sh`,
+then the ephemerides path should be automatically picked up from
+the $LALPULSAR_DATADIR environment variable.
+Alternatively, you can place a file
 `~/.pyfstat.conf` into your home directory which looks like
 
 ```
@@ -113,6 +137,9 @@ sun_ephem = '/home/<USER>/lalsuite-install/share/lalpulsar/sun00-19-DE421.dat.gz
 ```
 Paths set in this way will take precedence over the environment variable.
 
+Finally, you can manually specify ephemerides files when initialising
+each PyFstat search (as one of the arguments).
+
 ### Contributors
 
 * Greg Ashton

pyfstat/core.py

@@ -159,6 +159,8 @@ def predict_fstat(
     IFOs=None,
     assumeSqrtSX=None,
     tempory_filename="fs.tmp",
+    earth_ephem=None,
+    sun_ephem=None,
     **kwargs
 ):
     """ Wrapper to lalapps_PredictFstat
@@ -201,6 +203,11 @@ def predict_fstat(
     cl_pfs.append("--maxStartTime={}".format(int(maxStartTime)))
     cl_pfs.append("--outputFstat={}".format(tempory_filename))
 
+    if earth_ephem is not None:
+        cl_pfs.append("--ephemEarth='{}'".format(earth_ephem))
+    if sun_ephem is not None:
+        cl_pfs.append("--ephemSun='{}'".format(sun_ephem))
+
     cl_pfs = " ".join(cl_pfs)
     helper_functions.run_commandline(cl_pfs)
     d = read_par(filename=tempory_filename)
@@ -345,7 +352,8 @@ class BaseSearchClass(object):
             Paths of the two files containing positions of Earth and Sun,
             respectively at evenly spaced times, as passed to CreateFstatInput
 
-        Note: If not manually set, default values in ~/.pyfstat are used
+        Note: If not manually set, default values from get_ephemeris_files()
+              are used (looking in ~/.pyfstat or $LALPULSAR_DATADIR)
 
         """
 
@@ -353,8 +361,12 @@ class BaseSearchClass(object):
 
         if earth_ephem is None:
             self.earth_ephem = earth_ephem_default
+        else:
+            self.earth_ephem = earth_ephem
         if sun_ephem is None:
             self.sun_ephem = sun_ephem_default
+        else:
+            self.sun_ephem = sun_ephem
 
 
 class ComputeFstat(BaseSearchClass):

pyfstat/helper_functions.py

@@ -146,15 +146,24 @@ def get_ephemeris_files():
             earth_ephem = None
             sun_ephem = None
     elif env_var in list(os.environ.keys()):
-        earth_ephem = os.path.join(os.environ[env_var], "earth00-40-DE421.dat.gz")
-        sun_ephem = os.path.join(os.environ[env_var], "sun00-40-DE421.dat.gz")
+        ephem_version = "DE405"
+        earth_ephem = os.path.join(
+            os.environ[env_var], "earth00-40-{:s}.dat.gz".format(ephem_version)
+        )
+        sun_ephem = os.path.join(
+            os.environ[env_var], "sun00-40-{:s}.dat.gz".format(ephem_version)
+        )
         if not (os.path.isfile(earth_ephem) and os.path.isfile(sun_ephem)):
-            earth_ephem = os.path.join(os.environ[env_var], "earth00-19-DE421.dat.gz")
-            sun_ephem = os.path.join(os.environ[env_var], "sun00-19-DE421.dat.gz")
+            earth_ephem = os.path.join(
+                os.environ[env_var], "earth00-19-{:s}.dat.gz".format(ephem_version)
+            )
+            sun_ephem = os.path.join(
+                os.environ[env_var], "sun00-19-{:s}.dat.gz".format(ephem_version)
+            )
             if not (os.path.isfile(earth_ephem) and os.path.isfile(sun_ephem)):
                 logging.warning(
-                    "No [earth/sun]00-[19/40]-DE421 ephemerides "
-                    "found in the " + os.environ[env_var] + " directory. " + please
+                    "Default [earth/sun]00-[19/40]-" + ephem_version + " ephemerides "
+                    "not found in the " + os.environ[env_var] + " directory. " + please
                 )
                 earth_ephem = None
                 sun_ephem = None
@@ -368,3 +377,16 @@ def twoFDMoffThreshold(
         return twoFDMoffthreshold_below_threshold
     else:
         return 10 ** (prefactor * np.log10(twoFon - offset))
+
+
+def match_commandlines(cl1, cl2, be_strict_about_full_executable_path=False):
+    """ Check if two commandlines match element-by-element, regardless of order """
+    cl1s = cl1.split(" ")
+    cl2s = cl2.split(" ")
+    # first item will be the executable name
+    # by default be generous here and do not worry about full paths
+    if not be_strict_about_full_executable_path:
+        cl1s[0] = os.path.basename(cl1s[0])
+        cl2s[0] = os.path.basename(cl2s[0])
+    unmatched = np.setxor1d(cl1s, cl2s)
+    return len(unmatched) == 0

pyfstat/make_sfts.py

@@ -259,38 +259,65 @@ transientTau = {:10.0f}\n"""
     def check_cached_data_okay_to_use(self, cl_mfd):
         """ Check if cached data exists and, if it does, if it can be used """
 
-        getmtime = os.path.getmtime
+        need_new = "Will create new SFT file."
 
+        logging.info("Checking if cached data good to reuse...")
         if os.path.isfile(self.sftfilepath) is False:
-            logging.info("No SFT file matching {} found".format(self.sftfilepath))
-            return False
-        else:
-            logging.info("Matching SFT file found")
-
-        if getmtime(self.sftfilepath) < getmtime(self.config_file_name):
             logging.info(
-                (
-                    "The config file {} has been modified since the sft file {} "
-                    + "was created"
-                ).format(self.config_file_name, self.sftfilepath)
+                "No SFT file matching {} found. {}".format(self.sftfilepath, need_new)
             )
             return False
+        else:
+            logging.info("OK: Matching SFT file found.")
+
+        if "injectionSources" in cl_mfd:
+            if os.path.isfile(self.config_file_name):
+                if os.path.getmtime(self.sftfilepath) < os.path.getmtime(
+                    self.config_file_name
+                ):
+                    logging.info(
+                        (
+                            "The config file {} has been modified since the SFT file {} "
+                            + "was created. {}"
+                        ).format(self.config_file_name, self.sftfilepath, need_new)
+                    )
+                    return False
+                else:
+                    logging.info(
+                        "OK: The config file {} is older than the SFT file {}".format(
+                            self.config_file_name, self.sftfilepath
+                        )
+                    )
+                    # NOTE: at this point we assume it's safe to re-use, since
+                    # check_if_cff_file_needs_rewritting()
+                    # should have already been called before
+            else:
+                raise RuntimeError(
+                    "Commandline requires file '{}' but it is missing.".format(
+                        self.config_file_name
+                    )
+                )
 
-        logging.info(
-            "The config file {} is older than the sft file {}".format(
-                self.config_file_name, self.sftfilepath
-            )
-        )
-        logging.info("Checking contents of cff file")
+        logging.info("Checking new commandline against existing SFT header...")
         cl_dump = "lalapps_SFTdumpheader {} | head -n 20".format(self.sftfilepath)
         output = helper_functions.run_commandline(cl_dump)
-        found = [True for line in output.split("\n") if line[-len(cl_mfd) :] == cl_mfd]
-        if any(found):
-            logging.info("Contents matched, use old sft file")
-            return True
-        else:
-            logging.info("Contents unmatched, create new sft file")
+        header_lines_lalapps = [
+            line for line in output.split("\n") if "lalapps" in line
+        ]
+        if len(header_lines_lalapps) == 0:
+            logging.info(
+                "Could not obtain comparison commandline from old SFT header. "
+                + need_new
+            )
             return False
+        cl_old = header_lines_lalapps[0]
+        if not helper_functions.match_commandlines(cl_old, cl_mfd):
+            logging.info("Commandlines unmatched. " + need_new)
+            return False
+        else:
+            logging.info("OK: Commandline matched with old SFT header.")
+        logging.info("Looks like cached data matches current options, will re-use it!")
+        return True
 
     def check_if_cff_file_needs_rewritting(self, content):
         """ Check if the .cff file has changed
@@ -298,24 +325,27 @@ transientTau = {:10.0f}\n"""
         Returns True if the file should be overwritten - where possible avoid
         overwriting to allow cached data to be used
         """
+        logging.info("Checking if we can re-use injection config file...")
         if os.path.isfile(self.config_file_name) is False:
-            logging.info("No config file {} found".format(self.config_file_name))
+            logging.info("No config file {} found.".format(self.config_file_name))
             return True
         else:
-            logging.info("Config file {} already exists".format(self.config_file_name))
+            logging.info("Config file {} already exists.".format(self.config_file_name))
 
         with open(self.config_file_name, "r") as f:
             file_content = f.read()
             if file_content == content:
                 logging.info(
-                    "File contents match, no update of {} required".format(
+                    "File contents match, no update of {} required.".format(
                         self.config_file_name
                     )
                 )
                 return False
             else:
                 logging.info(
-                    "File contents unmatched, updating {}".format(self.config_file_name)
+                    "File contents unmatched, updating {}.".format(
+                        self.config_file_name
+                    )
                 )
                 return True
 
@@ -350,10 +380,16 @@ transientTau = {:10.0f}\n"""
         cl_mfd.append("--Tsft={}".format(self.Tsft))
         if self.h0 != 0:
             cl_mfd.append('--injectionSources="{}"'.format(self.config_file_name))
+        earth_ephem = getattr(self, "earth_ephem", None)
+        sun_ephem = getattr(self, "sun_ephem", None)
+        if earth_ephem is not None:
+            cl_mfd.append('--ephemEarth="{}"'.format(earth_ephem))
+        if sun_ephem is not None:
+            cl_mfd.append('--ephemSun="{}"'.format(sun_ephem))
 
         cl_mfd = " ".join(cl_mfd)
-
-        if self.check_cached_data_okay_to_use(cl_mfd) is False:
+        check_ok = self.check_cached_data_okay_to_use(cl_mfd)
+        if check_ok is False:
             helper_functions.run_commandline(cl_mfd)
 
     def predict_fstat(self):
@@ -371,6 +407,8 @@ transientTau = {:10.0f}\n"""
             self.detectors,
             self.sqrtSX,
             tempory_filename="{}.tmp".format(self.label),
+            earth_ephem=self.earth_ephem,
+            sun_ephem=self.sun_ephem,
         )  # detectors OR IFO?
         return twoF_expected
 
@@ -435,6 +473,7 @@ class GlitchWriter(Writer):
         see `lalapps_Makefakedata_v5 --help` for help with the other paramaters
         """
 
+        self.set_ephemeris_files()
         self.basic_setup()
         self.calculate_fmin_Band()
 
@@ -754,6 +793,12 @@ class FrequencyModulatedArtifactWriter(Writer):
         cl_mfd.append("--h0={}".format(h0))
         cl_mfd.append("--cosi={}".format(self.cosi))
         cl_mfd.append("--lineFeature=TRUE")
+        earth_ephem = getattr(self, "earth_ephem", None)
+        sun_ephem = getattr(self, "sun_ephem", None)
+        if earth_ephem is not None:
+            cl_mfd.append('--ephemEarth="{}"'.format(earth_ephem))
+        if sun_ephem is not None:
+            cl_mfd.append('--ephemSun="{}"'.format(sun_ephem))
         cl_mfd = " ".join(cl_mfd)
         helper_functions.run_commandline(cl_mfd, log_level=10)
 

tests.py

@@ -80,12 +80,16 @@ class Writer(Test):
         Writer = pyfstat.Writer(self.label, outdir=self.outdir, duration=3600)
         if os.path.isfile(Writer.sftfilepath):
             os.remove(Writer.sftfilepath)
+        # first run: make everything from scratch
         Writer.make_cff()
         Writer.run_makefakedata()
         time_first = os.path.getmtime(Writer.sftfilepath)
+        # second run: should re-use .cff and .sft
+        Writer.make_cff()
         Writer.run_makefakedata()
         time_second = os.path.getmtime(Writer.sftfilepath)
         self.assertTrue(time_first == time_second)
+        # third run: touch the .cff to force regeneration
         time.sleep(1)  # make sure timestamp is actually different!
         os.system("touch {}".format(Writer.config_file_name))
         Writer.run_makefakedata()