""" Classes needed for the known pulsar search pipeline. (C) 2006, 2015 Matthew Pitkin """ # make print statements python 3-proof from __future__ import print_function, division __author__ = 'Matthew Pitkin ' __date__ = '$Date$' __version__ = '$Revision$' import os import re from glue import pipeline import sys import ast import json import subprocess as sp import shutil import uuid from six.moves.configparser import RawConfigParser import six.moves.urllib.parse as urlparse from copy import deepcopy import numpy as np import pickle from scipy import optimize from collections import OrderedDict from six import string_types from lalapps import pulsarpputils as pppu # set some specific error codes and messages KNOPE_ERROR_GENERAL = -1 KNOPE_ERROR_NO_SEGMENTS = 2 KNOPE_ERROR_GENERAL_MSG = "Error... an error has occurred during DAG creation" KNOPE_ERROR_NO_SEGMENTS_MSG = "No required data segments were available to perform the analysis" # dictionary of error code/message pairs KNOPE_ERROR = {KNOPE_ERROR_GENERAL: KNOPE_ERROR_GENERAL_MSG, KNOPE_ERROR_NO_SEGMENTS: KNOPE_ERROR_NO_SEGMENTS_MSG} #set some specific warning codes KNOPE_WARNING_NO_SEGMENTS = 102 """ Class for setting up the DAG for the whole known pulsar search pipeline """ class knopeDAG(pipeline.CondorDAG): def __init__(self, cp, configfilename, pulsarlist=None): """ Initialise with ConfigParser cp object and the filename of the config file. If an error occurs the error_code variables will be set to -1. The value will stay as 0 on success. """ self.error_code = 0 # set error_code to KNOPE_ERROR_GENERAL following any failures self.warning_code = 0 self.config = cp if pulsarlist is not None: if isinstance(pulsarlist, list): self.pulsarlist = pulsarlist else: print("Error... 'pulsarlist' argument must be 'None' or a list.") self.error_code = KNOPE_ERROR_GENERAL return # if just re-doing post-processing try reading in previuos analysis pickle file self.postonly = self.get_config_option('analysis', 'postprocessing_only', cftype='boolean', default=False) prevdag = None if self.postonly: preprocessed_pickle = self.get_config_option('analysis', 'preprocessed_pickle_object') if preprocessed_pickle == None: print("Error... trying post-processing only, but no previous pickle file is given", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return else: try: fp = open(preprocessed_pickle, 'rb') prevdag = pickle.load(fp) fp.close() except: print("Error... trying post-processing only, but previous pickle file '%s' cannot be read in" % preprocessed_pickle, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return # Get run directory self.rundir = self.get_config_option('analysis', 'run_dir', cftype='dir', default=os.getcwd()) if self.error_code != 0: return # quit class if there's been an error # a list of allowed IFO names allowed_ifos = ['H1', 'H2', 'L1', 'G1', 'V1', 'T1', 'K1'] # Get the interferometers to analyse self.ifos = self.get_config_option('analysis', 'ifos', 'list') if self.error_code != 0: return # quit class if there's been an error # make sure ifo is an allowed value for ifo in self.ifos: if ifo not in allowed_ifos: print("Error... you have specified an unknown IFO '%s'" % ifo, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return if self.postonly: # check that this analysis uses the same, or a subset of, the previous IFOs for ifo in self.ifos: if ifo not in prevdag.ifos: print("Error... for 'post-processing-only' the current IFOs must be a subset of those in the previous run", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return # Get frequency factors (e.g. 1f and 2f) for analysis (default to twice the rotation frequency) self.freq_factors = self.get_config_option('analysis', 'freq_factors', cftype='list', default=[2.]) if len(self.freq_factors) > 2: # cannot have more than two values print("Warning... only up to two frequency factors can be given. Defaulting to [2.]") self.freq_factors = [2.] if len(self.freq_factors) == 2: # if there are two values they must be 1 and 2 if 1.0 not in self.freq_factors and 2.0 not in self.freq_factors: print("Warning... if giving two frequency factors they must be [1., 2.]. Defaulting to this") self.freq_factors = [1., 2.] for ff in self.freq_factors: if ff <= 0.: print("Warning... frequency factors cannot be negative. Defaulting to [2.]") self.freq_factors = [2.] if self.postonly: # check that this analysis uses the same, or a subset of, the previous frequency factors for ff in self.freq_factors: if ff not in prevdag.freq_factors: print("Error... for 'post-processing-only' the current frequency factors must be a subset of those in the previous run", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return # Get the base directory for the preprocessing analysis if not self.postonly: self.preprocessing_base_dir = self.get_config_option('analysis', 'preprocessing_base_dir', cftype='dict') if self.error_code != 0: return else: self.preprocessing_base_dir = prevdag.preprocessing_base_dir # try making directory/check if exists if not self.postonly: for ifo in self.ifos: self.mkdirs(self.preprocessing_base_dir[ifo]) if self.error_code != 0: return # Get the start and end time of the analysis # see if starttime and endtime are integers, or dictionaries (of integers or lists of integers) self.starttime = self.get_config_option('analysis', 'starttime') if self.error_code != 0: return self.starttime = ast.literal_eval(self.starttime) # check if int or dict if isinstance(self.starttime, int): # convert to dictionary stdict = {} for ifo in self.ifos: stdict[ifo] = [self.starttime] # convert to list self.starttime = stdict elif isinstance(self.starttime, dict): # check each detector has a start time for ifo in self.ifos: if ifo not in self.starttime: print("Error... 'starttime' either be a single 'int', or a dictionary containing all detectors with an integer or list of integers.", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return for stkey in dict(self.starttime): if isinstance(self.starttime[stkey], int): # convert to list self.starttime[stkey] = [self.starttime[stkey]] elif isinstance(self.starttime[stkey], list): # check all values are int if len([v for v in self.starttime[stkey] if isinstance(v, int)]) != len(self.starttime[stkey]): print("Error... 'starttime' either be a single 'int', or a dictionary containing all detectors with an integer or list of integers.", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return else: print("Error... 'starttime' either be a single 'int', or a dictionary containing all detectors with an integer or list of integers.", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return else: print("Error... 'starttime' either be a single 'int', or a dictionary containing all detectors with an integer or list of integers.", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return self.endtime = self.get_config_option('analysis', 'endtime') if self.error_code != 0: return self.endtime = ast.literal_eval(self.endtime) # check if int or dict if isinstance(self.endtime, int): # convert to dictionary etdict = {} for ifo in self.ifos: etdict[ifo] = [self.endtime] # convert to list self.endtime = etdict elif isinstance(self.endtime, dict): # check each detector has an end time for ifo in self.ifos: if ifo not in self.endtime: print("Error... 'endtime' either be a single 'int', or a dictionary containing all detectors with an integer or list of integers.", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return for etkey in dict(self.endtime): if isinstance(self.endtime[etkey], int): # convert to list self.endtime[etkey] = [self.endtime[etkey]] elif isinstance(self.endtime[etkey], list): # check all values are int if len([v for v in self.endtime[etkey] if isinstance(v, int)]) != len(self.endtime[etkey]): print("Error... 'endtime' either be a single 'int', or a dictionary containing all detectors with an integer or list of integers.", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return else: print("Error... 'endtime' either be a single 'int', or a dictionary containing all detectors with an integer or list of integers.", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return else: print("Error... 'endtime' either be a single 'int', or a dictionary containing all detectors with an integer or list of integers.", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return # check start time and end time lists are consistent self.ndatasets = {} # dictionary of number of seperate datasets (e.g. 2 different observing runs) for each detector for ifo in self.ifos: if len(self.starttime[ifo]) != len(self.endtime[ifo]): print("Error... 'starttime' and 'endtime' have an inconsistent number of entries.", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return else: self.ndatasets[ifo] = len(self.starttime[ifo]) # Get the pre-processing engine (heterodyne or SplInter - default to heterodyne) if not self.postonly: self.engine = self.get_config_option('analysis', 'preprocessing_engine', default='heterodyne') if self.error_code != 0: return if self.engine not in ['heterodyne', 'splinter']: print("Warning... 'preprocessing_engine' value '%s' not recognised. Defaulting to 'heterodyne'." % (self.engine)) self.engine = 'heterodyne' else: self.engine = prevdag.engine # Get the solar system ephemeris path try: defaultephempath = os.environ['LALPULSAR_DATADIR'] except: defaultephempath = None self.ephem_path = self.get_config_option('analysis', 'ephem_path', cftype='dir', default=defaultephempath) if self.error_code != 0: return # Get Condor accounting info (default to ligo.prod.o1.cw.targeted.bayesian) self.accounting_group = self.get_config_option('condor', 'accounting_group', default='ligo.prod.o1.cw.targeted.bayesian') if self.error_code != 0: return if 'cw.targeted.bayesian' not in self.accounting_group: print("Error... the 'accounting_group' should contain 'cw.targeted.bayesian'", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return if cp.has_option('condor', 'accounting_group_user'): self.accounting_group_user = cp.get('condor', 'accounting_group_user') if len(self.accounting_group_user) == 0: self.accounting_group_user = None # user not specified else: self.accounting_group_user = None # user not specified # Get the analysis run directory (where the DAG and sub files will be created and run from) (default to current working directory) self.run_dir = self.get_config_option('analysis', 'run_dir', default=os.getcwd()) self.mkdirs(self.run_dir) if self.error_code != 0: return # Get the analysis log directory self.log_dir = self.get_config_option('analysis', 'log_dir') self.mkdirs(self.log_dir) if self.error_code != 0: return uniqueid = str(uuid.uuid4().hex) daglog = self.get_config_option('analysis', 'dag_name', default='knope-'+uniqueid+'.log') if len(daglog) == 0: # if no dag_name was given, but dag_name was still in .ini file daglog = 'known_pulsar_pipeline-'+uniqueid+'.log' self.daglogfile = os.path.join(self.log_dir, daglog) # initialise DAG pipeline.CondorDAG.__init__(self, self.daglogfile) # set dag file dagname = self.get_config_option('analysis', 'dag_name', default='knope-'+uniqueid) if len(dagname) == 0: # if no dag_name was given, but dag_name was stillin .ini file dagname = 'known_pulsar_pipeline-'+uniqueid self.set_dag_file(os.path.join(self.run_dir, dagname)) # Check if running in autonomous mode self.autonomous = self.get_config_option('analysis', 'autonomous', cftype='boolean', default=False) if self.error_code != 0: return if self.autonomous: initialstart = self.get_config_option('analysis', 'autonomous_initial_start', cftype='int') if self.error_code != 0: return # convert to dictionary (to be consistent with start time) self.initial_start = {} for ifo in self.ifos: self.initial_start[ifo] = [initialstart] # convert to list to be consistent with self.starttimes self.initial_start = self.initial_start else: self.initial_start = self.starttime # Get pulsars to analyse (in the future this should be able to get pulsars from a database based on # certain parameters) self.pulsar_param_dir = self.get_config_option('analysis', 'pulsar_param_dir') if not os.path.isdir(self.pulsar_param_dir) and not os.path.isfile(self.pulsar_param_dir): print("Error... pulsar parameter file/directory '%s' does not exist!" % self.pulsar_param_dir, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return if self.postonly: if prevdag.pulsar_param_dir != self.pulsar_param_dir: print("Error... for 'post-processing-only' the pulsar parameter directory must be that same as in the previous run", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return self.unmodified_pulsars = [] # a list of pulsar .par files that have not been modified since the last run self.modified_pulsars = [] # a list of pulsar .par files that are new, or have been modified since the last run self.modification_files = [] # a list of dictionaries (one for each pulsar) containing the modification time file and modification time # Go through pulsar directory and check for any modified/new files if os.path.isdir(self.pulsar_param_dir): self.param_files = [os.path.join(self.pulsar_param_dir, pf) for pf in os.listdir(self.pulsar_param_dir) if '.par' in pf and os.path.isfile(os.path.join(self.pulsar_param_dir, pf)) and '.mod_' not in pf] if len(self.param_files) == 0: print("Error... no pulsar parameter files found in '%s'" % self.pulsar_param_dir, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return self.param_files.sort() # sort the files into alphabetical order elif os.path.isfile(self.pulsar_param_dir): # if a single file convert into list self.param_files = [self.pulsar_param_dir] else: print("Error... pulsar parameter file or directory '%s' does not exist" % self.pulsar_param_dir, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return self.analysed_pulsars = {} # dictionary keyed by PSRJ and containing the .par file path for all analysed pulsars self.skipped_pulsars = {} # dictionary keyed on .par files for all skipped pulsars (including reason for exclusion) if not self.postonly: for par in self.param_files: # check that there is a PSRJ name in the .par file - if not then skip this pulsar # Get par file data try: psr = pppu.psr_par(par) except: print("Could not read in parameter file '%s'. Skipping this pulsar." % par) self.skipped_pulsars[par] = "Could not read in parameter file" continue # check that there is a PSRJ name in the .par file (this is required) if 'PSRJ' not in psr.__dict__: print("Could not read 'PSRJ' value from '%s'. Skipping this pulsar." % par) self.skipped_pulsars[par] = "Could not read 'PSRJ' value" continue # check if we're just using selected input pulsars if pulsarlist is not None: if psr['PSRJ'] not in pulsarlist: continue # check that (if the pulsar is in a binary) it is an allowed binary type (given what's coded in lalpulsar) if 'BINARY' in psr.__dict__: bintype = psr['BINARY'] if bintype not in ['BT', 'BT1P', 'BT2P', 'BTX', 'ELL1', 'DD', 'DDS', 'MSS', 'T2']: print("Binary type '%s' in '%s' is not recognised. Skipping this pulsar." % (bintype, par)) self.skipped_pulsars[par] = "Binary type '%s' is currenlty not recognised" % bintype continue # check that .par file contains ephemeris information and units (if not present defaults will be used - see get_ephemeris) if psr['EPHEM'] != None: if psr['EPHEM'] not in ['DE200', 'DE405', 'DE414', 'DE421', 'DE430', 'DE435', 'DE436']: print("Unregconised ephemeris '%s' in '%s'. Skipping this source" % (psr['EPHEM'], par)) self.skipped_pulsars[par] = "Unregconised ephemeris '%s'" % psr['EPHEM'] continue if psr['UNITS'] != None: if psr['UNITS'] not in ['TCB', 'TDB']: print("Unregconised time units '%s' in '%s'. Skipping this source" % (psr['UNITS'], par)) self.skipped_pulsars[par] = "Unregconised ephemeris '%s'" % psr['UNITS'] continue self.analysed_pulsars[psr['PSRJ']] = par # create parfile modification time file modtimefile = os.path.join(os.path.dirname(par), '.mod_'+os.path.basename(par)) # Check if modification time file exists if not os.path.isfile(modtimefile): # add par file to modified pulsars list self.modified_pulsars.append(par) # create modification time file modtime = os.stat(par).st_mtime self.modification_files.append({'file': modtimefile, 'time': str(modtime)}) else: # check whether par file modification time is consistent with value in modtimefile parmodtime = str(os.stat(par).st_mtime) fm = open(modtimefile, 'r') try: oldmodtime = fm.readline().strip() except: print("Warning... could not read modification time from '%s'. Assuming file is modified" % modtimefile) oldmodtime = -1.23456789 fm.close() if parmodtime == oldmodtime: # file is unmodified self.unmodified_pulsars.append(par) else: # file is modified self.modification_files.append({'file': modtimefile, 'time': parmodtime}) # update the time in the .mod file self.modified_pulsars.append(par) if pulsarlist is not None: if len(self.analysed_pulsars) == 0: print("Could not find any of the listed pulsars '[%s]' in the .par file directory '%s'." % (', '.join(pulsarlist), self.pulsar_param_dir)) self.error_code = KNOPE_ERROR_GENERAL return self.segment_file_update = [] # to contain a list of pairs of segment files - the former to be concatenated into the latter if required if not self.postonly: # Setup pre-processing jobs (datafind, segment finding and heterodyne/splinter running) for each pulsar self.setup_preprocessing() if self.error_code != 0: return # Check whether only the preprocessing (heterodyne or splinter) jobs are required self.preprocessing_only = self.get_config_option('analysis', 'preprocessing_only', cftype='boolean', default=False) if self.error_code != 0: return else: # use information from previous run self.preprocessing_only = False self.remove_job = None if pulsarlist is None: self.unmodified_pulsars = prevdag.unmodified_pulsars self.modified_pulsars = prevdag.modified_pulsars self.analysed_pulsars = prevdag.analysed_pulsars self.skipped_pulsars = prevdag.skipped_pulsars self.processed_files = prevdag.processed_files else: # if analysing only selected pulsars (from pulsarlist) just get information on them self.processed_files = {} for psrl in pulsarlist: if psrl not in prevdag.analysed_pulsars: print("Error... specified pulsar '%s' could not be found in previous run pickle file '%s'." % (psrl, preprocessed_pickle)) self.error_code = KNOPE_ERROR_GENERAL return if psrl in prevdag.unmodified_pulsars: self.unmodified_pulsars[psrl] = prevdag.unmodified_pulsars[psrl] if psrl in prevdag.modified_pulsars: self.modified_pulsars[psrl] = prevdag.modified_pulsars[psrl] self.analysed_pulsars[psrl] = prevdag.analysed_pulsars[psrl] self.processed_files[psrl] = prevdag.processed_files[psrl] if self.preprocessing_only: # end class initialisation # output par file modification times (do this now, so that modification files are overwritten if the script had failed earlier) for pitem in self.modification_files: fm = open(pitem['file'], 'w') fm.write(pitem['time']) fm.close() # output amended segment list files for sfs in self.segment_file_update: p = sp.Popen("cat " + sfs[0] + " >> " + sfs[1], shell=True) p.communicate() if p.returncode != 0: print("Warning... could not append segments to previous segments file. No log of previous segments will be available.") return # setup parameter estimation self.setup_parameter_estimation() if self.error_code != 0: return # setup post-processing page creation self.setup_results_pages() ### FINAL CLOSE OUT ITEMS # output par file modification times (do this now, so that modification files are overwritten if the script had failed earlier) for pitem in self.modification_files: fm = open(pitem['file'], 'w') fm.write(pitem['time']) fm.close() # output ammended segment list files for sfs in self.segment_file_update: p = sp.Popen("cat " + sfs[0] + " >> " + sfs[1], shell=True) p.communicate() if p.returncode != 0: print("Warning... could not append segments to previous segments file. No log of previous segments will be available.") # output dictionaries of analysed pulsars and skipped pulsars to JSON files fpa = open(os.path.join(self.run_dir, 'analysed_pulsars.txt'), 'w') json.dump(self.analysed_pulsars, fpa, indent=2) fpa.close() fps = open(os.path.join(self.run_dir, 'skipped_pulsars.txt'), 'w') json.dump(self.skipped_pulsars, fps, indent=2) fps.close() # email notification that the analysis has finished if required email = self.get_config_option('analysis', 'email', allownone=True) if email != None: if '@' not in email: print("Warning... email address '%s' is invalid. No notification will be sent." % email) else: import smtplib import socket # try setting email server try: # set up sender try: HOST = socket.getfqdn() USER = os.environ['USER'] FROM = USER+'@'+HOST except: # default from to 'matthew.pitkin@ligo.org' FROM = 'matthew.pitkin@ligo.org' subject = "lalapps_knope: successful setup" messagetxt = "Hi User,\n\nYour analysis using configuration file '%s' has successfully setup the analysis. Once complete the results will be found at %s.\n\nRegards\n\nlalapps_knope\n" % (configfilename, self.results_url) emailtemplate = "From: {0}\nTo: {1}\nSubject: {2}\n\n{3}" message = emailtemplate.format(FROM, email, subject, messagetxt) server = smtplib.SMTP('localhost') server.sendmail(FROM, email, message) server.quit() except: print("Warning... could not send notification email.") def setup_results_pages(self): """ Setup the results webpage creation """ # get directory for output results self.results_basedir = self.get_config_option('results_page', 'web_dir') self.results_pulsar_dir = {} # dictionary of results directories for individual pulsars self.results_pulsar_ini = {} # dictionary of results configuration files for individual pulsars self.results_pulsar_url = {} # dictionary of results page URLs self.mkdirs(self.results_basedir) if self.error_code != 0: return # get URL for output results self.results_url = self.get_config_option('results_page', 'base_url') # run individual pulsar results page creation self.results_exec = self.get_config_option('results_page', 'results_exec', default='/usr/bin/lalapps_knope_result_page') self.results_universe = self.get_config_option('results_page', 'universe', default='local') # check file exists and is executable if not os.path.isfile(self.results_exec) or not os.access(self.results_exec, os.X_OK): print("Warning... 'results_exec' in '[results_page]' does not exist or is not an executable. Try finding code in path.") resultexec = self.find_exec_file('lalapps_knope_result_page') if resultexec == None: print("Error... could not find 'lalapps_knope_result_page' in 'PATH'", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return else: self.results_exec = resultexec self.collate_exec = self.get_config_option('results_page', 'collate_exec', default='/usr/bin/lalapps_knope_collate_results') # check file exists and is executable if not os.path.isfile(self.collate_exec) or not os.access(self.collate_exec, os.X_OK): print("Warning... 'collate_exec' in '[results_page]' does not exist or is not an executable. Try finding code in path.") collateexec = self.find_exec_file('lalapps_knope_collate_results') if collateexec == None: print("Error... could not find 'lalapps_knope_collate_results' in 'PATH'", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return else: self.collate_exec = collateexec # check if running on an injection self.injection = self.get_config_option('analysis', 'injections', cftype='boolean', default=False) self.use_gw_phase = self.get_config_option('pe', 'use_gw_phase', cftype='boolean', default=False) if self.error_code != 0: return # check upper limit credible interval to use self.upper_limit = self.get_config_option('results_page', 'upper_limit', cftype='int', default=95) # check whether to show joint posterior plot for all parameters self.show_all_posteriors = self.get_config_option('results_page', 'show_all_posteriors', cftype='boolean', default=False) # check whether to subtract injected/heterodyned values from phase parameters for plots self.subtract_truths = self.get_config_option('results_page', 'subtract_truths', cftype='boolean', default=False) # check whether to plot priors on 1D posteriors plots self.show_priors = self.get_config_option('results_page', 'show_priors', cftype='boolean', default=False) # check whether to copy "all" files used to create results (fine heterodyne files, par file, prior file, posterior files) # into the results page directory self.copy_all_files = self.get_config_option('results_page', 'copy_all_files', cftype='boolean', default=False) cpjob = None if self.copy_all_files: # create job for copying files cpjob = copyJob(accgroup=self.accounting_group, accuser=self.accounting_group_user, logdir=self.log_dir, rundir=self.run_dir) # create parameter estimation job resultpagejob = resultpageJob(self.results_exec, univ=self.results_universe, accgroup=self.accounting_group, accuser=self.accounting_group_user, logdir=self.log_dir, rundir=self.run_dir) collatejob = collateJob(self.collate_exec, univ=self.results_universe, accgroup=self.accounting_group, accuser=self.accounting_group_user, logdir=self.log_dir, rundir=self.run_dir) # create config file for collating results into a results table cpc = RawConfigParser() # create config parser to output .ini file # create configuration .ini file cinifile = os.path.join(self.results_basedir, 'collate.ini') # create sections cpc.add_section('output') cpc.add_section('input') cpc.add_section('general') cpc.set('output', 'path', self.results_basedir) # set output directory cpc.set('input', 'path', self.results_basedir) # set directory containing individual results directories # get sort type and direction sorttype = self.get_config_option('results_page', 'sort_value', default='name') # default sorting on name cpc.set('general', 'sort_value', sorttype) sortdirection = self.get_config_option('results_page', 'sort_direction', default='ascending') # default sorting in ascending order cpc.set('general', 'sort_direction', sortdirection) if self.pe_coherent_only: cpc.set('general', 'detectors', ['Joint']) elif self.pe_incoherent_only: cpc.set('general', 'detectors', self.ifos) else: cdets = deepcopy(self.ifos) cdets.append('Joint') cpc.set('general', 'detectors', cdets) # get pulsar parameters to output paramout = self.get_config_option('results_page', 'parameters', cftype='list', default=['f0']) cpc.set('general', 'parameters', paramout) # get results to output resout = self.get_config_option('results_page', 'results', cftype='list', default=['h0ul']) cpc.set('general', 'results', resout) # write and set ini file try: fp = open(cinifile, 'w') cpc.write(fp) fp.close() except: print("Error... could not write configuration file '%s' for results collation page" % cinifile, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return # loop through pulsars for pname in self.analysed_pulsars: # create results directory for pulsar self.results_pulsar_dir[pname] = os.path.join(self.results_basedir, pname) self.results_pulsar_url[pname] = urlparse.urljoin(self.results_url, pname) self.mkdirs(self.results_pulsar_dir[pname]) if self.error_code != 0: return copydir = None if self.copy_all_files: # create directory for file copies copydir = os.path.join(self.results_pulsar_dir[pname], 'files') self.mkdirs(copydir) if self.error_code != 0: return # copy par file shutil.copy(self.analysed_pulsars[pname], os.path.join(copydir, pname+'.par')) # copy prior file shutil.copy(self.pe_prior_files[pname], copydir) # copy correlation coefficient file if pname in self.pe_cor_files: shutil.copy(self.pe_cor_files[pname], copydir) # if in autonomous mode, and a previous JSON results file already exists, make a copy of it jsonfile = os.path.join(self.results_pulsar_dir[pname], pname+'.json') if self.autonomous: if os.path.isfile(jsonfile): # append starttime (which will be the end time of the previous results) timestamp to JSON file try: shutil.copyfile(jsonfile, jsonfile + '_%d' % list(self.starttime.values())[0]) except: print("Warning... could not copy previous results JSON file '%s'. Previous results may get overwritten." % jsonfile, file=sys.stderr) # try getting some pulsar information (dist, p1_I and assoc) from the ATNF catalogue for use in lalapps_knope_result_page # NOTE: this is mainly required because on the ARCCA cluster the nodes cannot access the internet pinfo = pppu.get_atnf_info(pname) if pinfo is not None: dist, p1_I, assoc, _ = pinfo # unpack values psrinfo = {} psrinfo['Pulsar data'] = {} psrinfo['Pulsar data']['DIST'] = dist psrinfo['Pulsar data']['P1_I'] = p1_I psrinfo['Pulsar data']['ASSOC'] = assoc try: fp = open(jsonfile, 'w') json.dump(psrinfo, fp, indent=2) fp.close() except: print("Warning... could not write out ATNF catalogue information to JSON file '%s'." % jsonfile, file=sys.stderr) cp = RawConfigParser() # create config parser to output .ini file # create configuration .ini file inifile = os.path.join(self.results_pulsar_dir[pname], pname+'.ini') # add sections cp.add_section('general') cp.add_section('parameter_estimation') cp.add_section('data') cp.add_section('output') cp.add_section('plotting') cp.set('output', 'path', self.results_pulsar_dir[pname]) # set output directory cp.set('output', 'indexpage', os.path.relpath(self.results_basedir, self.results_pulsar_dir[pname])) cp.set('general', 'parfile', self.analysed_pulsars[pname]) # set the pulsar parameter file cp.set('general', 'detectors', self.ifos) # set the detectors cp.set('general', 'upper_limit', self.upper_limit) # set the upper limit credible interval if self.pe_coherent_only: cp.set('general', 'joint_only', True) # only output the joint multi-detector analysis else: cp.set('general', 'joint_only', False) if self.pe_incoherent_only: cp.set('general', 'with_joint', False) # only output individual detector analyses else: cp.set('general', 'with_joint', True) # include joint multi-detector analysis if self.pe_num_background > 0: cp.set('general', 'with_background', True) # include background analysis else: cp.set('general', 'with_background', False) # no background analysis present if self.injection: cp.set('general', 'injection', True) # set if an injection or not else: cp.set('general', 'injection', False) if self.use_gw_phase: cp.set('general', 'use_gw_phase', True) # use GW initial phase (rather than rotational phase) e.g. for hardware injections else: cp.set('general', 'use_gw_phase', False) cp.set('general', 'harmonics', self.freq_factors) # set frequency harmonics in analysis cp.set('general', 'model_type', self.pe_model_type) # set 'waveform' or 'source' model type cp.set('general', 'biaxial', self.pe_biaxial) # set if using a biaxial source model if self.show_priors: cp.set('general', 'priorfile', self.pe_prior_files[pname]) # set the prior file # get posterior files (and background directories) posteriorsfiles = {} backgrounddir = {} copydirpos = None if copydir is not None: copydirpos = os.path.join(copydir, 'posteriors') self.mkdirs(copydirpos) for i, comb in enumerate(self.pe_combinations): dets = comb['detectors'] detprefix = comb['prefix'] if len(dets) == 1: det = dets[0] else: det = 'Joint' # use 'Joint' as the term for a multi-detector analysis posteriorsfiles[det] = os.path.join(self.pe_posterior_basedir, pname) posteriorsfiles[det] = os.path.join(posteriorsfiles[det], detprefix) if copydirpos is not None: self.mkdirs(os.path.join(copydirpos, detprefix)) if self.pe_num_background > 0: backgrounddir[det] = os.path.join(self.pe_posterior_background_basedir, pname) backgrounddir[det] = os.path.join(backgrounddir[det], detprefix) dirpostfix = '' if len(self.freq_factors) > 1: # coherent multi-frequency analysis dirpostfix = 'multiharmonic' else: if not self.freq_factors[0]%1.: # for integers just output directory as e.g. 2f dirpostfix = '%df' % int(self.freq_factors[0]) else: dirpostfix = '%.2ff' % self.freq_factors[0] posteriorsfiles[det] = os.path.join(posteriorsfiles[det], dirpostfix) posteriorsfiles[det] = os.path.join(posteriorsfiles[det], 'posterior_samples_%s.hdf' % pname) if copydirpos is not None: copydirposp = os.path.join(copydirpos, detprefix) copydirposp = os.path.join(copydirposp, dirpostfix) self.mkdirs(copydirposp) cpnode = copyNode(cpjob) cpnode.set_source(posteriorsfiles[det]) cpnode.set_destination(copydirposp) for n2pnode in self.pe_nest2pos_nodes[pname]: cpnode.add_parent(n2pnode) self.add_node(cpnode) # if in autonomous mode copy previous posterior files if self.autonomous: if os.path.isfile(posteriorsfiles[det]): try: # copy to file with the start time (i.e. the end time of the previous analysis for which the posterior file belongs) appended shutil.copyfile(posteriorsfiles[det], posteriorsfiles[det].strip('.hdf') + '_%d.hdf' % list(self.starttime.values())[0]) except: print("Warning... could not create copy of current posterior samples file '%s'. This will get overwritten on next autonomous run." % posteriorsfiles[det], file=sys.stderr) if self.pe_num_background > 0: backgrounddir[det] = os.path.join(backgrounddir[det], dirpostfix) cp.set('parameter_estimation', 'posteriors', posteriorsfiles) if self.pe_num_background > 0: cp.set('parameter_estimation', 'background', backgrounddir) # copy fine heterodyned/spectrally interpolated files if required copydirhet = None if copydir is not None: copydirhet = os.path.join(copydir, 'data') self.mkdirs(copydirhet) datafiles = {} for ifo in self.ifos: if copydir is not None: copydirhet = os.path.join(copydirhet, ifo) self.mkdirs(copydirhet) filelist = [] for ff in self.freq_factors: filelist.append(self.processed_files[pname][ifo][ff][-1]) # copy fine heterodyned/spectrally interpolated files if copydir is not None: if not ff%1.: # for integers just output director as e.g. 2f ffdir = os.path.join(copydirhet, '%df' % int(ff)) else: # for non-integers use 2 d.p. for dir name ffdir = os.path.join(copydirhet, '%.3ff' % int(ff)) self.mkdirs(ffdir) cpnode = copyNode(cpjob) cpnode.set_source(self.processed_files[pname][ifo][ff][-1]) cpnode.set_destination(ffdir) if self.engine == 'heterodyne' and not self.postonly: if pname in self.concat_nodes and self.autonomous: cpnode.add_parent(self.concat_nodes[pname][ifo][ff]) else: cpnode.add_parent(self.fine_heterodyne_nodes[pname][ifo][ff]) elif self.engine == 'splinter' and not self.postonly: cpnode.add_parent(self.splinter_nodes_modified[ifo][ff]) self.add_node(cpnode) if len(self.freq_factors) == 1: filelist = filelist[0] # don't need list if only one value datafiles[ifo] = filelist if copydir is not None: copydirhet = os.path.join(copydir, 'data') # reset to base path cp.set('data', 'files', datafiles) cp.set('plotting', 'all_posteriors', self.show_all_posteriors) cp.set('plotting', 'subtract_truths', self.subtract_truths) # output configuration file try: fp = open(inifile, 'w') cp.write(fp) fp.close() except: print("Error... could not write configuration file '%s' for results page" % inifile, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return # create dag node resultsnode = resultpageNode(resultpagejob) resultsnode.set_config(inifile) # add parent lalinference_nest2pos job for n2pnode in self.pe_nest2pos_nodes[pname]: resultsnode.add_parent(n2pnode) if self.pe_num_background > 0: for n2pnode in self.pe_nest2pos_background_nodes[pname]: resultsnode.add_parent(n2pnode) self.add_node(resultsnode) # add as parent to the collation node collatenode = collateNode(collatejob) # create a collate node for each pulsar, so that the table gets regenerated as each new reults comes in collatenode.set_config(cinifile) collatenode.add_parent(resultsnode) self.add_node(collatenode) def setup_preprocessing(self): """ Setup the preprocessing analysis: data finding, segment finding and heterodyne/splinter data processing """ # set the data find job and nodes self.setup_datafind() if self.error_code != 0: return # loop through the pulsars and setup required jobs for each if self.engine == 'heterodyne': self.setup_heterodyne() if self.error_code != 0: return if self.engine == 'splinter': self.setup_splinter() if self.error_code != 0: return # create jobs to concatenate output fine heterodyned/Splinter files if required self.remove_job = None self.concatenate_files() if self.error_code != 0: return def setup_parameter_estimation(self): """ Setup parameter estimation jobs/nodes for signal and background analyses """ # get executable self.pe_exec = self.get_config_option('pe', 'pe_exec', default='lalapps_pulsar_parameter_estimation_nested') if self.error_code != 0: return # check file exists and is executable if not os.path.isfile(self.pe_exec) or not os.access(self.pe_exec, os.X_OK): print("Warning... 'pe_exec' in '[pe]' does not exist or is not an executable. Try finding code in path.") peexec = self.find_exec_file('lalapps_pulsar_parameter_estimation_nested') if peexec == None: print("Error... could not find 'lalapps_pulsar_parameter_estimation_nested' in 'PATH'", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return else: self.pe_exec = peexec # Get Condor universe (default to vanilla) self.pe_universe = self.get_config_option('pe', 'universe', default='vanilla') if self.error_code != 0: return self.pe_nest2pos_nodes = {} # condor nodes for the lalinference_nest2pos jobs (needed to use as parents for results processing jobs) self.pe_nest2pos_background_nodes = {} # nodes for background analysis # see whether to run just as independent detectors, or independently AND coherently over all detectors if len(self.ifos) == 1: self.pe_incoherent_only = True # set to True for one detector else: self.pe_incoherent_only = self.get_config_option('analysis', 'incoherent_only', cftype='boolean', default=False) # see whether to run only the coherent multidetector analysis analysis self.pe_coherent_only = False if len(self.ifos) > 1: self.pe_coherent_only = self.get_config_option('analysis', 'coherent_only', cftype='boolean', default=False) # get the number of background analyses to perform self.pe_num_background = self.get_config_option('analysis', 'num_background', cftype='int', default=0) if self.pe_num_background < 0: print("Warning... 'num_background' is a negative value. Defaulting to zero background runs") self.pe_num_background = 0 # get the PE output directory self.pe_output_basedir = self.get_config_option('pe', 'pe_output_dir', cftype='dir') # Make directory self.mkdirs(self.pe_output_basedir) if self.error_code != 0: return # set background run directories if required self.pe_output_background_basedir = self.get_config_option('pe', 'pe_output_dir_background', cftype='dir', allownone=True) if self.pe_num_background != 0: if self.pe_output_background_basedir == None: print("Error... no background analysis directory has been set", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL; return else: self.mkdirs(self.pe_output_background_basedir) if self.error_code != 0: return # get some general PE parameters self.pe_nlive = self.get_config_option('pe', 'n_live', cftype='int', default=2048) # number of live points self.pe_nruns = self.get_config_option('pe', 'n_runs', cftype='int', default=1) # number of parallel runs self.pe_tolerance = self.get_config_option('pe', 'tolerance', cftype='float', default=0.1) # nested sampling stopping criterion self.pe_random_seed = self.get_config_option('pe', 'random_seed', cftype='int', allownone=True) # random number generator seed self.pe_nmcmc = self.get_config_option('pe', 'n_mcmc', cftype='int', allownone=True) # number of MCMC steps for each nested sample self.pe_nmcmc_initial = self.get_config_option('pe', 'n_mcmc_initial', cftype='int', default=500) self.pe_non_gr = self.get_config_option('pe', 'non_gr', cftype='boolean', default=False) # use non-GR parameterisation (default to False) # check if only using a section of data self.pe_starttime = self.get_config_option('pe', 'starttime', cftype='float', allownone=True) self.pe_endtime = self.get_config_option('pe', 'endtime', cftype='float', allownone=True) self.pe_truncate_time = self.get_config_option('pe', 'truncate_time', cftype='float', allownone=True) self.pe_truncate_samples = self.get_config_option('pe', 'truncate_samples', cftype='int', allownone=True) self.pe_truncate_fraction = self.get_config_option('pe', 'truncate_fraction', cftype='float', allownone=True) # parameters for background runs self.pe_nruns_background = self.get_config_option('pe', 'n_runs_background', cftype='int', default=1) self.pe_nlive_background = self.get_config_option('pe', 'n_live_background', cftype='int', default=1024) # parameters for ROQ self.pe_roq = self.get_config_option('pe', 'use_roq', cftype='boolean', default=False) # check if using Reduced Order Quadrature (ROQ) self.pe_roq_ntraining = self.get_config_option('pe', 'roq_ntraining', cftype='int', default=2500) # number of training waveforms for reduced basis generation self.pe_roq_tolerance = self.get_config_option('pe', 'roq_tolerance', cftype='float', default=5e-12) # mis-match tolerance when producing reduced basis self.pe_roq_uniform = self.get_config_option('pe', 'roq_uniform', cftype='boolean', default=False) # check if setting uniform distributions for sprinkling traning waveform parameters self.pe_roq_chunkmax = self.get_config_option('pe', 'roq_chunkmax', cftype='int', default=1440) # maximum data chunk length for creating ROQ # FIXME: Currently this won't run with non-GR parameters, so output a warning and default back to GR! if self.pe_non_gr: print("Warning... currently this will not run with non-GR parameters. Reverting to GR-mode.") self.pe_non_gr = False # if searching at both the rotation frequency and twice rotation frequency set which parameterisation to use self.pe_model_type = self.get_config_option('pe', 'model_type', default='waveform') if self.pe_model_type not in ['waveform', 'source']: print("Warning... the given 'model_type' '%s' is not allowed. Defaulting to 'waveform'" % self.pe_model_type) self.pe_model_type = 'waveform' self.pe_biaxial = False # check whether using a biaxial signal if len(self.freq_factors) == 2 and self.pe_non_gr == False: self.pe_biaxial = self.get_config_option('pe', 'biaxial', cftype='boolean', default=False) # use a biaxial signal model # check whether using the Student's t-likelihood or Gaussian likelihood self.pe_gaussian_like = self.get_config_option('pe', 'gaussian_like', cftype='boolean', default=False) if self.engine == 'splinter': # always use the Gaussian likelihood for the SplInter processed data self.pe_gaussian_like = True # check if there is a pre-made prior file to use for all pulsars self.pe_prior_options = self.get_config_option('pe', 'prior_options', cftype='dict') self.pe_premade_prior_file = self.get_config_option('pe', 'premade_prior_file', allownone=True) if self.pe_premade_prior_file is not None: if not os.path.isfile(self.pe_premade_prior_file): print("Error... pre-made prior file '{}' does not exist!".format(self.pe_premade_prior_file), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return self.pe_derive_amplitude_prior = self.get_config_option('pe', 'derive_amplitude_prior', cftype='boolean', default=False) if self.pe_derive_amplitude_prior: # get JSON file containing any previous amplitude upper limits for pulsars self.pe_amplitude_prior_file = self.get_config_option('pe', 'amplitude_prior_file', allownone=True) # get JSON file containing any previous posterior files from which to use amplitude vs cosiota samples as GMM prior self.pe_previous_posterior_files = None self.pe_previous_posteriors_file = self.get_config_option('pe', 'previous_posteriors_file', allownone=True) self.pe_prior_info = None self.pe_prior_asds = {} # get file, or dictionary of amplitude spectral density files (e.g. from a previous run) to derive amplitude priors try: self.pe_amplitude_prior_asds = ast.literal_eval(self.get_config_option('pe', 'amplitude_prior_asds', allownone=True)) except: self.pe_amplitude_prior_asds = self.get_config_option('pe', 'amplitude_prior_asds', allownone=True) try: self.pe_amplitude_prior_obstimes = ast.literal_eval(self.get_config_option('pe', 'amplitude_prior_obstimes', allownone=True)) except: self.pe_amplitude_prior_obstimes = self.get_config_option('pe', 'amplitude_prior_obstimes', allownone=True) self.pe_amplitude_prior_type = self.get_config_option('pe', 'amplitude_prior_type', default='fermidirac') # check if using par files errors in parameter estimation self.pe_use_parameter_errors = self.get_config_option('pe', 'use_parameter_errors', cftype='boolean', default=False) # check for executable for merging nested sample files/converting them to posteriors self.pe_n2p_exec = self.get_config_option('pe', 'n2p_exec', default='lalinference_nest2pos') if self.error_code != 0: return # check file exists and is executable if not os.path.isfile(self.pe_n2p_exec) or not os.access(self.pe_n2p_exec, os.X_OK): print("Warning... 'pe_n2p_exec' in '[pe]' does not exist or is not an executable. Try finding code in path.") pen2pexec = self.find_exec_file('lalinference_nest2pos') if pen2pexec == None: print("Error... could not find 'lalinference_nest2pos' in 'PATH'", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return else: self.pe_n2p_exec = pen2pexec self.pe_posterior_basedir = self.get_config_option('pe', 'n2p_output_dir') if self.pe_posterior_basedir == None: print("Error... no 'n2p_output_dir' specified in '[pe]' giving path for posterior sample outputs", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return self.mkdirs(self.pe_posterior_basedir) if self.error_code != 0: return if self.pe_num_background > 0: self.pe_posterior_background_basedir = self.get_config_option('pe', 'n2p_output_dir_background') if self.pe_posterior_background_basedir == None: print("Error... no 'n2p_output_dir_background' specified in '[pe]' giving path for posterior sample outputs", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return self.mkdirs(self.pe_posterior_background_basedir) if self.error_code != 0: return # check whether to remove all the nested sample files after the run self.pe_clean_nest_samples = self.get_config_option('pe', 'clean_nest_samples', cftype='boolean', default=False) # check if there are memory requirements for the PE jobs self.pe_request_memory = self.get_config_option('pe', 'pe_request_memory', cftype='int', allownone=True) # create parameter estimation job pejob = ppeJob(self.pe_exec, univ=self.pe_universe, accgroup=self.accounting_group, accuser=self.accounting_group_user, logdir=self.log_dir, rundir=self.run_dir, requestmemory=self.pe_request_memory) # create job to merge nested samples and convert to posterior samples n2pjob = nest2posJob(self.pe_n2p_exec, univ='local', accgroup=self.accounting_group, accuser=self.accounting_group_user, logdir=self.log_dir, rundir=self.run_dir) # create job for moving SNR files mvjob = moveJob(accgroup=self.accounting_group, accuser=self.accounting_group_user, logdir=self.log_dir, rundir=self.run_dir) # create job for removing nested samples (use previous remove job if existing) if self.remove_job == None: rmjob = removeJob(accgroup=self.accounting_group, accuser=self.accounting_group_user, logdir=self.log_dir, rundir=self.run_dir) else: rmjob = self.remove_job # work out combinations of parameter estimation jobs to run (FIXME: if non-GR mode was able to be enabled then another # combination would include those below, but in both GR and non-GR flavours) # NOTE: if running with two frequency factors (which have been fixed, so they can only be 1f and 2f), this just runs # for a multi-harmonic search at both frequencies. If you want to run at a single frequency then that should be set up # independently. self.pe_combinations = [] for ifo in self.ifos: if not self.pe_coherent_only: self.pe_combinations.append({'detectors': [ifo], 'prefix': ifo}) # all individual detector runs if not self.pe_incoherent_only: self.pe_combinations.append({'detectors': self.ifos, 'prefix': "".join(self.ifos)}) # add joint multi-detector run self.pe_prior_files = {} # dictionary to contain prior files for each pulsar self.pe_cor_files = {} # dictionary to contain correlation coefficient matrices for pulsars as required # dictionary of nest2pos nodes n2pnodes = {} # loop over the total number of jobs (single signal job and background jobs) njobs = self.pe_num_background + 1 for j in range(njobs): # loop through each pulsar that has been analysed for pname in self.analysed_pulsars: psr = pppu.psr_par(self.analysed_pulsars[pname]) # directories for output nested sample files if j == 0: # first iteration of outer loop is the signal job psrdir = os.path.join(self.pe_output_basedir, pname) else: # subsequent iterations are for background jobs psrdir = os.path.join(self.pe_output_background_basedir, pname) # directories for output posterior files if j == 0: # first iteration of outer loop is the signal job psrpostdir = os.path.join(self.pe_posterior_basedir, pname) else: # subsequent iterations are for background jobs psrpostdir = os.path.join(self.pe_posterior_background_basedir, pname) if self.pe_roq: nroqruns = 1 # add one run for the ROQ weights calculation else: nroqruns = 0 # create directory for that pulsar self.mkdirs(psrdir) if self.error_code != 0: return self.mkdirs(psrpostdir) if self.error_code != 0: return n2pnodes[pname] = [] # list of nodes for lalinference_nest2pos jobs for a given pulsar for comb in self.pe_combinations: dets = comb['detectors'] detprefix = comb['prefix'] # set up directories for detectors detdir = os.path.join(psrdir, detprefix) self.mkdirs(detdir) if self.error_code != 0: return detpostdir = os.path.join(psrpostdir, detprefix) self.mkdirs(detpostdir) if self.error_code != 0: return # make directories for frequency factors if len(self.freq_factors) > 1: # coherent multi-frequency analysis ffdir = os.path.join(detdir, 'multiharmonic') ffpostdir = os.path.join(detpostdir, 'multiharmonic') else: if not self.freq_factors[0]%1.: # for integers just output directory as e.g. 2f ffdir = os.path.join(detdir, '%df' % int(self.freq_factors[0])) ffpostdir = os.path.join(detpostdir, '%df' % int(self.freq_factors[0])) else: ffdir = os.path.join(detdir, '%.2ff' % self.freq_factors[0]) ffpostdir = os.path.join(detpostdir, '%.2ff' % self.freq_factors[0]) # add integer numbering to directories if running background analysis jobs if j > 0: ffdir = os.path.join(ffdir, '%05d' % (j-1)) # start directory names as 00000, 00001, 00002, etc ffpostdir = os.path.join(ffpostdir, '%05d' % (j-1)) self.mkdirs(ffdir) if self.error_code != 0: return self.mkdirs(ffpostdir) if self.error_code != 0: return randomiseseed = '' if j == 0: # create prior file for analysis (use this same file for all background runs) priorfile = self.create_prior_file(psr, psrdir, dets, self.freq_factors, ffdir) if pname not in self.pe_prior_files: self.pe_prior_files[pname] = priorfile # set prior file (just use first one as they should be the same for each combination of detectors) nruns = self.pe_nruns nlive = self.pe_nlive else: nruns = self.pe_nruns_background nlive = self.pe_nlive_background # set seed for randomising data (this needs to be the same over each nruns) randomiseseed = ''.join([str(f) for f in np.random.randint(1, 10, size=15).tolist()]) # set output ROQ weights file if self.pe_roq: roqweightsfile = os.path.join(ffdir, 'roqweights.bin') nestfiles = [] # list of nested sample file names penodes = [] roqinputnode = None # setup job(s) i = counter = 0 while counter < nruns+nroqruns: # loop over the required number of runs penode = ppeNode(pejob, psrname=pname) if self.pe_random_seed is not None: penode.set_randomseed(self.pe_random_seed) # set seed for RNG penode.set_detectors(','.join(dets)) # add detectors penode.set_par_file(self.analysed_pulsars[pname]) # add parameter file if pname in self.pe_cor_files: penode.set_cor_file(self.pe_cor_files[pname]) # set parameter correlation coefficient file penode.set_prior_file(priorfile) # set prior file penode.set_harmonics(','.join([str(ffv) for ffv in self.freq_factors])) # set frequency harmonics penode.set_Nlive(nlive) # set number of live points if self.pe_nmcmc is not None: penode.set_Nmcmc(self.pe_nmcmc) # set number of MCMC iterations for choosing new points penode.set_Nmcmcinitial(self.pe_nmcmc_initial) # set initial number of MCMC interations for reampling the prior penode.set_tolerance(self.pe_tolerance) # set tolerance for ending nested sampling if self.pe_starttime is not None: penode.set_start_time(self.pe_starttime) # set start time of data to use if self.pe_endtime is not None: penode.set_end_time(self.pe_endtime) # set end time of data to use if self.pe_starttime is None and self.pe_endtime is None: if self.pe_truncate_time is not None: penode.set_truncate_time(self.pe_truncate_time) elif self.pe_truncate_samples is not None: penode.set_truncate_samples(self.pe_truncate_samples) elif self.pe_truncate_fraction is not None: penode.set_truncate_fraction(self.pe_truncate_fraction) # set the output nested samples file nestfiles.append(os.path.join(ffdir, 'nested_samples_%s_%05d.hdf' % (pname, i))) penode.set_outfile(nestfiles[i]) if self.pe_roq: penode.set_roq() penode.set_roq_chunkmax(str(self.pe_roq_chunkmax)) if counter == 0: # first time round just output weights penode.set_roq_ntraining(str(self.pe_roq_ntraining)) penode.set_roq_tolerance(str(self.pe_roq_tolerance)) penode.set_roq_outputweights(roqweightsfile) if self.pe_roq_uniform: penode.set_roq_uniform() else: # use pre-created weights file penode.set_roq_inputweights(roqweightsfile) # for background runs set the randomise seed if j > 0: penode.set_randomise(randomiseseed) # add input files inputfiles = [] for det in dets: for ff in self.freq_factors: inputfiles.append(self.processed_files[pname][det][ff][-1]) penode.set_input_files(','.join(inputfiles)) if self.pe_gaussian_like or self.engine == 'splinter': # use Gaussian likelihood penode.set_gaussian_like() if (len(self.freq_factors) == 2 and 1. in self.freq_factors and 2. in self.freq_factors): # set whether using source model if self.pe_model_type == 'source': penode.set_source_model() # set whether using a biaxial signal if self.pe_biaxial: penode.set_biaxial() # set Earth, Sun and time ephemeris files earthfile, sunfile, timefile = self.get_ephemeris(psr) penode.set_ephem_earth(earthfile) penode.set_ephem_sun(sunfile) penode.set_ephem_time(timefile) # add parents (unless just doing post-processing) if not self.postonly: for ff in self.freq_factors: for det in dets: # set parents of node if pname in self.concat_nodes: penode.add_parent(self.concat_nodes[pname][det][ff]) else: if self.engine == 'heterodyne': try: # fine heterodyne nodes might not exist if (in automated mode) no new segments were found, but previous data was available penode.add_parent(self.fine_heterodyne_nodes[pname][det][ff]) except: pass if self.engine == 'splinter': try: # splinter nodes might not exist if (in automated mode) no new segments were found, but previous data was available if pname in self.splinter_modified_pars: penode.add_parent(self.splinter_nodes_modified[det][ff]) elif pname in self.splinter_unmodified_pars: penode.add_parent(self.splinter_nodes_unmodified[det][ff]) except: pass # if using ROQ add first PE node as parent to the rest if self.pe_roq: if roqinputnode is not None: # add first penode (which generates the ROQ interpolant) as a parent to subsequent nodes penode.add_parent(roqinputnode) if counter == 0: # get first penode (generating and outputting the ROQ interpolant) to add as parents to subsequent nodes roqinputnode = penode self.add_node(penode) counter = counter+1 # increment "counter", but not "i" del nestfiles[-1] # remove last element of nestfiles continue # add node to dag self.add_node(penode) penodes.append(penode) # move SNR files into posterior directory mvnode = moveNode(mvjob) snrsourcefile = os.path.splitext(nestfiles[i])[0]+'_SNR' # source SNR file snrdestfile = os.path.join(ffpostdir, 'SNR_%05d.txt' % i) # destination SNR file in posterior directory mvnode.set_source(snrsourcefile) mvnode.set_destination(snrdestfile) mvnode.add_parent(penode) self.add_node(mvnode) counter = counter+1 i = i+1 # add lalinference_nest2pos node to combine outputs/convert to posterior samples n2pnode = nest2posNode(n2pjob) postfile = os.path.join(ffpostdir, 'posterior_samples_%s.hdf' % pname) n2pnode.set_outfile(postfile) # output posterior file n2pnode.set_nest_files(nestfiles) # nested sample files n2pnodes[pname].append(n2pnode) # add parents for pn in penodes: n2pnode.add_parent(pn) self.add_node(n2pnode) # add node # check whether to remove the nested sample files if self.pe_clean_nest_samples: rmnode = removeNode(rmjob) # add name of header file rmnode.set_files(nestfiles) rmnode.add_parent(n2pnode) self.add_node(rmnode) if j == 0: # add main n2p nodes for pname in self.analysed_pulsars: self.pe_nest2pos_nodes[pname] = n2pnodes[pname] else: # add background n2p nodes for pname in self.analysed_pulsars: self.pe_nest2pos_background_nodes[pname] = n2pnodes[pname] def create_prior_file(self, psr, psrdir, detectors, freqfactors, outputpath, scalefactor=25.): """ Create the prior file to use for a particular job defined by a set of detectors, or single detector, and a set of frequency factors, or a single frequency factor. If creating an prior limit based on a set of given amplitude spectral densities (by calculating an estimate of the 95% UL they would produce) then it will additionally be scaled by a factor of `scalefactor`. Return the full output file and the create prior node """ # create the output file pname = psr['PSRJ'] outfile = os.path.join(outputpath, '%s.prior' % pname) # if using a pre-made prior file then just create a symbolic link to that file into outputpath if self.pe_premade_prior_file is not None: try: os.symlink(self.pe_premade_prior_file, outfile) except: print("Error... could not create symbolic link to prior file '%s'" % self.pe_premade_prior_file, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile # check if requiring to add parameters with errors in the .par file to the prior options prior_options = {} if self.pe_prior_options is not None: prior_options = deepcopy(self.pe_prior_options) if self.pe_use_parameter_errors: # create and output a covariance matrix (in the pulsar directory) if it does not already exist corfile = os.path.join(psrdir, '%s.cor' % pname) fp = open(corfile, 'w') fp.write('\t') # indent the list of parameters erritems = [] # list of values with errors ignore_pars = ["DM", "START", "FINISH", "NTOA", "TRES", "TZRMJD", "TZRFRQ", "TZRSITE", "NITS", "ELAT", "ELONG"] # keys to ignore from par file for paritem in pppu.float_keys: if paritem in ignore_pars: continue if psr['%s_ERR' % paritem] != None and psr['%s_FIT' % paritem] != None: # get values with a given error (suffixed with _ERR) if psr['%s_FIT' % paritem] == 1: # set Gaussian prior with mean being the parameter value and sigma being the error if paritem in ['RA_RAD', 'DEC_RAD']: itemname = paritem.replace('_RAD', '') else: itemname = paritem prior_options[itemname] = {'priortype': 'gaussian', 'ranges': [psr[paritem], psr['%s_ERR' % paritem]]} fp.write(itemname+' ') erritems.append(itemname) if len(erritems) > 0: fp.write('\n') for i, ei in enumerate(erritems): # have all values uncorrelated except w0 and T0 or wdot and Pb, which should be highly correlated for small eccentricity binaries fp.write(ei+'\t') for j in range(i+1): if i == j: fp.write('1 ') # diagonal values of correlation coefficient matrix else: # check if an eccentricity of between 0 and 0.001 ecc = 0.0 if psr['E'] != None: ecc = psr['E'] elif psr['ECC'] != None: ecc = psr['ECC'] if ecc >= 0. and ecc < 0.001: if ((ei == 'T0' and erritems[j] == 'OM') or (ei == 'OM' and erritems[j] == 'T0')) and ('T0' in erritems and 'OM' in erritems): fp.write('0.9999 ') # set parameters to be highly correlated (although not fully due to issues with fully correlated parameters) elif ((ei == 'PB' and erritems[j] == 'OMDOT') or (ei == 'OMDOT' and erritems[j] == 'PB')) and ('PB' in erritems and 'OMDOT' in erritems): fp.write('0.9999 ') # set parameters to be highly correlated (although not fully due to issues with fully correlated parameters) else: fp.write('0 ') else: fp.write('0 ') # set everything else to be uncorrelated fp.write('\n') fp.close() if len(erritems) > 0: self.pe_cor_files[pname] = corfile else: # no error value were found so remove corfile os.remove(corfile) # check if deriving amplitude priors if self.pe_derive_amplitude_prior: # check if there is a file containing a previous posterior to use for amplitude vs cosiota prior posteriorfile = None if self.pe_previous_posterior_files is None: if self.pe_previous_posteriors_file is not None: try: fp = open(self.pe_previous_posteriors_file, 'r') self.pe_previous_posterior_files = json.load(fp) fp.close() posteriorfile = self.pe_previous_posterior_files[pname] except: print("Error... could not open file '%s' listing previous posterior files." % self.pe_previous_posteriors_file, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile else: # check if current pulsar has a previous posterior sample file to use as prior if psr in self.pe_previous_posterior_files: posteriorfile = self.pe_previous_posterior_files[pname] # open output prior file try: fp = open(outfile, 'w') except: print("Error... could not open prior file '%s'" % outfile, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile # write out any priors that have been given for prioritem in prior_options: if 'priortype' not in prior_options[prioritem]: print("Error... no 'priortype' given for parameter '%s'" % prioritem, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile ptype = prior_options[prioritem]['priortype'] if ptype != 'gmm': if 'ranges' not in prior_options[prioritem]: print("Error... no 'ranges' given for parameter '%s'" % prioritem, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile rangevals = prior_options[prioritem]['ranges'] if len(rangevals) != 2: print("Error... 'ranges' for parameter '%s' must be a list or tuple with two entries" % prioritem, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile # ignore cosiota if using previous posterior file as prior if posteriorfile is not None and prioritem.upper() == 'COSIOTA': continue else: fp.write('%s\t%s\t%.16le\t%.16le\n' % (prioritem, ptype, rangevals[0], rangevals[1])) else: # check if item is prior for multiple parameters npars = len(prioritem.split(':')) # output if a Gaussian Mixture Model prior is set if 'nmodes' not in prior_options[prioritem]: print("Error... no 'nmodes' given for parameter '{}'".format(prioritem), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile nmodes = prior_options[prioritem]['nmodes'] if 'means' not in prior_options[prioritem]: print("Error... no 'means' given for parameter '{}'".format(prioritem), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile means = prior_options[prioritem]['means'] if len(means) != nmodes: print("Error... number of mean values must be equal to the number of modes", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile for mean in means: if len(mean) != npars: print("Error... number of mean values must be equal to the number of parameters", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile if 'covs' not in prior_options[prioritem]: print("Error... no 'covs' given for parameter '{}'".format(prioritem), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile covs = prior_options[prioritem]['covs'] if len(means) != nmodes: print("Error... number of covariance matrices values must be equal to the number of modes", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile for cov in covs: npcov = np.array(cov) if npcov.shape[0] != npcov.shape[1] and npcov.shape[1] != npars: print("Error... number of covariance matrices rows/columns must be equal to the number of parameters", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile if 'weights' not in prior_options[prioritem]: print("Error... no 'weights' given for parameter '{}'".format(prioritem), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile weights = prior_options[prioritem]['weights'] if len(weights) != nmodes: print("Error... number of weights must be equal to the number of modes", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile if 'ranges' in prior_options[prioritem]: ranges = prior_options[prioritem]['ranges'] if len(ranges) != npars: print("Error... number of ranges must be equal to the number of parameters", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile for rangevals in ranges: if len(rangevals) != 2: print("Error... ranges must have two values", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile else: ranges = None fp.write('{}\tgmm\t'.format(prioritem)) fp.write('{}\t'.format(nmodes)) fp.write('{}\t'.format(re.sub(r'\s+', '', str(means)))) fp.write('{}\t'.format(re.sub(r'\s+', '', str(covs)))) fp.write('{}'.format(re.sub(r'\s+', '', str(weights)))) if ranges is not None: for rangevals in ranges: fp.write('\t') fp.write('{}'.format(re.sub(r'\s+', '', str(rangevals)))) fp.write('\n') # set the required amplitude priors or parameters to estimate from previous posterior samples (and limits) requls = {} # dictionary to contain the required upper limits gmmpars = OrderedDict() if self.pe_model_type == 'waveform': if 2. in self.freq_factors: requls['C22'] = None gmmpars['C22'] = [0., np.inf] # limits on C22 prior if 1. in self.freq_factors: requls['C21'] = None gmmpars['C21'] = [0., np.inf] # limits on C21 prior elif self.pe_model_type == 'source': if len(self.freq_factors) == 1: requls['H0'] = None gmmpars['H0'] = [0., np.inf] # limits on H0 prior if len(self.freq_factors) == 2: if 1. in self.freq_factors and 2. in self.freq_factors: requls['I21'] = None requls['I31'] = None gmmpars['I21'] = [0., np.inf] gmmpars['I31'] = [0., np.inf] gmmpars['COSIOTA'] = [-1., 1.] # limits on COSIOTA prior (required in all cases) if len(requls) == 0: print("Error... unknown frequency factors or model type in configuration file.", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile # try and get the file containing previous upper limits if os.path.isfile(self.pe_amplitude_prior_file): # check file can be read if self.pe_prior_info is None: try: fpp = open(self.pe_amplitude_prior_file, 'r') self.pe_prior_info = json.load(fpp) # should be JSON file fpp.close() except: print("Error... could not parse prior file '%s'." % self.pe_amplitude_prior_file, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile # see if pulsar is in prior file if pname in self.pe_prior_info: uls = self.pe_prior_info[pname] for ult in requls: if ult == 'C22': if 'C22UL' not in uls and 'H0UL' in uls: # use 'H0' value for 'C22' if present requls['C22'] = uls['H0UL'] else: if ult+'UL' in uls: requls[ult] = uls[ult+'UL'] # if there are some required amplitude limits that have not been obtained try and get amplitude spectral densities freq = psr['F0'] if None in requls.values() and freq > 0.0 and posteriorfile is None: if self.pe_amplitude_prior_asds is not None and self.pe_amplitude_prior_obstimes is not None: asdfiles = self.pe_amplitude_prior_asds obstimes = self.pe_amplitude_prior_obstimes if not isinstance(asdfiles, dict): # if a single file is given convert into dictionary asdfilestmp = {} obstimestmp = {} asdfilestmp['det'] = asdfiles obstimestmp['det'] = float(obstimes) asdfiles = asdfilestmp obstimes = obstimestmp asdlist = [] for dk in asdfiles: # read in all the ASD files if dk not in obstimes: print("Error... no corresponding observation times for detector '%s'" % dk, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile else: if not isinstance(obstimes[dk], float) and not isinstance(obstimes[dk], int): print("Error... observation time must be a float or int.", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile if dk not in self.pe_prior_asds: if not os.path.isfile(asdfiles[dk]): print("Error... ASD file '%s' does not exist." % asdfiles[dk], file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile else: try: self.pe_prior_asds[dk] = np.loadtxt(asdfiles[dk], comments=['%', '#']) except: print("Error... could not load file '%s'." % asdfiles[dk], file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile asd = self.pe_prior_asds[dk] asdv = [] # empty array if 1. in self.freq_factors and (asd[0,0] <= freq and asd[-1,0] >= freq): # add ASD at 1f idxf = (np.abs(asd[:,0]-freq)).argmin() # get value nearest required frequency asdv.append(asd[idxf,1]) if 2. in self.freq_factors and (asd[0,0] <= 2.*freq and asd[-1,0] >= 2.*freq): idxf = (np.abs(asd[:,0]-2.*freq)).argmin() # get value nearest required frequency asdv.append(asd[idxf,1]) if len(asdv) > 0: asdlist.append(np.array(asdv)**2/(obstimes[dk]*86400.)) else: print("Error... frequency range in ASD file does not span pulsar frequency.", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile # get upper limit spectrum (harmonic mean of all the weighted spectra) mspec = np.zeros(len(self.freq_factors)) for asdv in asdlist: # interpolate frequencies mspec = mspec + (1./asdv) mspec = np.sqrt(1./mspec) # final weighted spectrum ulspec = 10.8*mspec # scaled to given "averaged" 95% upper limit estimate # set upper limits for creating priors if self.pe_model_type == 'waveform': if 1. in self.freq_factors: if requls['C21'] == None: requls['C21'] = ulspec[self.freq_factors.index(1.0)]*scalefactor if 2. in self.freq_factors: if requls['C22'] == None: requls['C22'] = ulspec[self.freq_factors.index(2.0)]*scalefactor if self.pe_model_type == 'source': if len(self.freq_factors) == 1: if requls['H0'] == None: requls['H0'] = ulspec[0]*scalefactor else: if 1. in self.freq_factors and 2. in self.freq_factors: # set both I21 and I31 to use the maximum of the 1f and 2f estimate if requls['I21'] == None: requls['I21'] = np.max(ulspec)*scalefactor if requls['I31'] == None: requls['I31'] = np.max(ulspec)*scalefactor # get amplitude prior type if self.pe_amplitude_prior_type not in ['fermidirac', 'uniform'] and posteriorfile is None: print("Error... prior type must be 'fermidirac' or 'uniform'", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile if posteriorfile is None: # go through required upper limits and output a Fermi-Dirac prior that also has a 95% limit at that value for ult in requls: if requls[ult] == None: print("Error... a required upper limit for '%s' is not available." % ult, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile else: if self.pe_amplitude_prior_type == 'fermidirac': try: b, a = self.fermidirac_rsigma(requls[ult]) except: print("Error... problem deriving the Fermi-Dirac prior for '%s'." % ult, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile else: a = 0. # uniform prior bound at 0 b = requls[ult]/0.95 # stretch limit to ~100% bound fp.write('%s\t%s\t%.16le\t%.16le\n' % (ult, self.pe_amplitude_prior_type, a, b)) else: # try and fit Gaussian Mixture Model to required amplitude parameters means, covs, weights, _ = self.gmm_prior(posteriorfile, gmmpars, taper='elliptical', decaywidth=1.) if self.error_code == -1: print("Error... could not set GMM prior using previous posterior samples.", file=sys.stderr) return outfile # output GMM prior parssep = ':'.join(gmmpars.keys()) fp.write("%s\tgmm\t%d\t" % (parssep, len(means))) # write out means meanstr = ','.join(['['+','.join([str(v) for v in vs.tolist()])+']' for vs in means]) fp.write("[%s]\t" % meanstr) # write out covariance matrices covstr = ','.join(['['+','.join(['['+','.join([str(ca) for ca in c])+']' for c in cs.tolist()])+']' for cs in covs]) fp.write("[%s]\t" % covstr) # write out weights fp.write("[%s]\t" % ','.join([str(w) for w in weights])) # write out limits for each parameter in turn for gp in gmmpars: fp.write("[%s]\t" % ','.join([str(lim) for lim in gmmpars[gp]])) fp.write("\n") fp.close() else: # make prior file from parameters try: fp = open(outfile, 'w') except: print("Error... could not write prior file '%s'" % outfile, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile for prioritem in prior_options: ptype = prior_options[prioritem]['priortype'] if ptype != 'gmm': rangevals = prior_options[prioritem]['ranges'] if len(rangevals) != 2: print("Error... the ranges in the prior for '%s' are not set properly" % prioritem, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile fp.write("%s\t%s\t%.16e\t%.16e\n" % (prioritem, ptype, rangevals[0], rangevals[1])) else: # check if item is prior for multiple parameters npars = len(prioritem.split(':')) # output if a Gaussian Mixture Model prior is set if 'nmodes' not in prior_options[prioritem]: print("Error... no 'nmodes' given for parameter '{}'".format(prioritem), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile nmodes = prior_options[prioritem]['nmodes'] if 'means' not in prior_options[prioritem]: print("Error... no 'means' given for parameter '{}'".format(prioritem), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile means = prior_options[prioritem]['means'] if len(means) != nmodes: print("Error... number of mean values must be equal to the number of modes", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile for mean in means: if len(mean) != npars: print("Error... number of mean values must be equal to the number of parameters", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile if 'covs' not in prior_options[prioritem]: print("Error... no 'covs' given for parameter '{}'".format(prioritem), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile covs = prior_options[prioritem]['covs'] if len(means) != nmodes: print("Error... number of covariance matrices values must be equal to the number of modes", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile for cov in covs: npcov = np.array(cov) if npcov.shape[0] != npcov.shape[1] and npcov.shape[1] != npars: print("Error... number of covariance matrices rows/columns must be equal to the number of parameters", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile if 'weights' not in prior_options[prioritem]: print("Error... no 'weights' given for parameter '{}'".format(prioritem), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile weights = prior_options[prioritem]['weights'] if len(weights) != nmodes: print("Error... number of weights must be equal to the number of modes", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile if 'ranges' in prior_options[prioritems]: ranges = prior_options[prioritem]['ranges'] if len(ranges) != npars: print("Error... number of ranges must be equal to the number of parameters", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile for rangevals in ranges: if len(rangevals) != 2: print("Error... ranges must have two values", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return outfile else: ranges = None fp.write('{}\tgmm\t'.format(prioritem)) fp.write('{}\t'.format(nmodes)) fp.write('{}\t'.format(re.sub(r'\s+', '', str(means)))) fp.write('{}\t'.format(re.sub(r'\s+', '', str(covs)))) fp.write('{}'.format(re.sub(r'\s+', '', str(weights)))) if ranges is not None: for rangevals in ranges: fp.write('\t') fp.write('{}'.format(re.sub(r'\s+', '', str(rangevals)))) fp.write('\n') fp.close() return outfile def fermidirac_rsigma(self, ul, mufrac=0.4, cdf=0.95): """ Calculate the r and sigma parameter of the Fermi-Dirac distribution to be used. Based on the definition of the distribution given in https://www.authorea.com/users/50521/articles/65214/_show_article the distribution will be defined by a mu parameter at which the distribution has 50% of it's maximum probability, and mufrac which is the fraction of mu defining the range from which the distribution falls from 97.5% of the maximum down to 2.5%. Using an upper limit defining a given cdf of the distribution the parameters r and sigma will be returned. """ Z = 7.33 # factor that defined the 97.5% -> 2.5% probability attenuation band around mu r = 0.5*Z/mufrac # set r # using the Fermi-Dirac CDF to find sigma given a distribution where the cdf value is found at ul solution = optimize.root(lambda s: cdf*np.log(1.+np.exp(r))-np.log(1.+np.exp(-r))-(ul/s)-np.log(1.+np.exp((ul/s)-r)), ul) sigma = solution.x[0] return r, sigma def gmm_prior(self, prevpostfile, pardict, ncomps=20, taper=None, decaywidth=5.): """ Create an ND Gaussian Mixture Model for use as a prior. This will use the BayesianGaussianMixture Model from scikit-learn, which fits a Dirichlet Process Gaussian Mixture Model to the input data infering the number of components required. The input to this should be a previously calculated posterior sample file, or numpy array of samples. If a files is given then the parameters given as keys in the `pardict` ordered dictionary will be extracted. For each parameter name key in the `pardict` ordered there should be pairs of hard upper an lower limits of the particular parameters. If any of these are not +/-infinity then the samples will be duplicated and reflected around that limit. This is to avoid edge effects for the inferred Gaussian distributions. `ncomps` sets the hyperparameter used in the Dirichlet process related to the number of Gaussian components. `taper` sets whether or not to taper-off any reflected samples, and how that tapering happens. Tapering can use: a 'gaussian' taper, where the half-width of the Gaussian is set by the range of the samples multiplied by `decaywidth`; a 'triangular' taper, which falls from one to zero over the range of the samples; an 'exponential' taper, where the decay constant is defined by 'decaywidth' multiplied by the range of the samples; or, an 'elliptical' taper, where the axis of the ellipse is set by 'decaywidth' multiplied by the range of the samples. The default is that no tapering is applied, and it should be noted that tapering can still leave artifacts in the final GMM. The means, covariance matrices and weights of the Gaussian components will be returned, along with the full set of points (including reflected points) used for the estimation. An example of using this would be for "H0" versus "COSIOTA", in which case the `pardict` might be: >> pardict = OrderedDict() >> pardict['H0'] = [0., np.inf] >> pardict['COSIOTA'] = [-1., 1.] """ try: from sklearn import mixture except: print('Error... could not import scikit-learn.', file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return None, None, None, None means = None covs = None weights = None if not isinstance(pardict, OrderedDict): print('Error... Input must be an ordered dictionary', file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return means, covs, weights, None npars = len(pardict) allsamples = [] # get samples try: if not isinstance(prevpostfile, (np.ndarray, np.generic)): # get samples from posterior sample hdf5 file if not os.path.isfile(prevpostfile): print("Error... previous posterior sample file '%s' does not exist" % prevpostfile, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return means, covs, weights, None possamps, _, _ = pppu.pulsar_nest_to_posterior(prevpostfile) for par in pardict: allsamples.append(possamps[par.upper()].samples) else: # get samples fron numpy array if prevpostfile.shape[1] == npars: for i in range(npars): allsamples.append(prevpostfile[:,i].reshape(len(prevpostfile), 1)) else: print('Error... input numpy array does not contain correct number of parameters', file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return means, covs, weights, None except: print("Error... could not extract posterior samples from file or numpy array", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return means, covs, weights, None # reflect and duplicate samples if required for all parameters (for each reflection add to ncomp) allsamplesnp = np.copy(allsamples).squeeze().T for i, p in enumerate(pardict.keys()): refsamples = None for lim in pardict[p]: if np.isfinite(lim): maxp = np.max(allsamples[i]) minp = np.min(allsamples[i]) sigmap = decaywidth*(maxp-minp) dist = lim - allsamplesnp[:,i] refidxs = np.ones(len(allsamplesnp[:,i]), dtype=bool) # reflect about this limit (with given tapering if required) if taper is not None: if lim > maxp: deltav = allsamplesnp[:,i]+2.*dist-maxp elif lim < minp: deltav = minp-(allsamplesnp[:,i]+2.*dist) else: print("Warning... limit is inside the extent of the samples", file=sys.stderr) continue probkeep = np.ones(len(allsamplesnp[:,i])) if taper == 'gaussian': probkeep = np.exp(-0.5*(deltav)**2/sigmap**2) elif taper == 'triangular': probkeep = 1.-(deltav)/(maxp-minp) elif taper == 'exponential': probkeep = np.exp(-(deltav)/sigmap) elif taper == 'elliptical': probkeep = np.zeros(len(allsamplesnp[:,i])) probkeep[deltav < sigmap] = np.sqrt(1.-(deltav[deltav < sigmap]/sigmap)**2) else: print("Warning... unknown tapering has been set, so none will be applied", file=sys.stderr) refidxs = (np.random.rand(len(allsamplesnp[:,i])) < probkeep).flatten() thesesamples = allsamplesnp[refidxs,:] thesesamples[:,i] += 2.*dist[refidxs] if refsamples is None: refsamples = np.copy(thesesamples) else: refsamples = np.concatenate((refsamples, thesesamples)) # stack samples if refsamples is not None: allsamplesnp = np.concatenate((allsamplesnp, refsamples)) # scale parameters to avoid dynamic range issues parscales = np.std(allsamplesnp, axis=0) scalesmat = np.identity(npars)*parscales scaledsamples = allsamplesnp/parscales # run DPGMM (tolerance and max_iter are increased over the default values # to aid convergence, although "unconverged" initialisations are generally # still fine to use) dpgmm = mixture.BayesianGaussianMixture(n_components=ncomps, covariance_type='full', tol=5e-2, max_iter=500).fit(scaledsamples) # predict the GMM components in which the samples live parpred = dpgmm.predict(scaledsamples) # get the means, covariances and weights of the GMM components in which actually contain predicted samples means = [] covs = [] weights = [] for i, (mean, covar, weight) in enumerate(zip(dpgmm.means_, dpgmm.covariances_, dpgmm.weights_)): if np.any(parpred == i): # check if any samples predicted to be in component i # check that mode is within 3.5 sigma of limits otherwise don't include it outlims = 0 for (mus, sigs, lowlim, highlim) in zip(mean*parscales, parscales*np.sqrt(np.diag(covar)), [pardict[p][0] for p in pardict], [pardict[p][1] for p in pardict]): if mus < lowlim - 3.*sigs or mus > highlim + 3.*sigs: outlims += 1 if outlims == 2: continue # rescale to get back to true values means.append(mean*parscales) covs.append(np.dot(scalesmat, np.dot(covar, scalesmat))) weights.append(weight) if len(means) == 0: print("Error... no GMM components returned", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return means, covs, weights, allsamplesnp def setup_heterodyne(self): """ Setup the coarse and fine heterodyne jobs/nodes. """ # get executable self.heterodyne_exec = self.get_config_option('heterodyne', 'heterodyne_exec', default='lalapps_heterodyne_pulsar') if self.error_code != 0: return # check file exists and is executable if not os.path.isfile(self.heterodyne_exec) or not os.access(self.heterodyne_exec, os.X_OK): print("Warning... 'heterodyne_exec' in '[heterodyne]' does not exist or is not an executable. Try finding code in path.") hetexec = self.find_exec_file('lalapps_heterodyne_pulsar') if hetexec == None: print("Error... could not find 'lalapps_heterodyne_pulsar' in 'PATH'", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return else: self.heterodyne_exec = hetexec # Get Condor universe (default to vanilla) self.heterodyne_universe = self.get_config_option('heterodyne', 'universe', default='vanilla') if self.error_code != 0: return # Get some general setup data self.coarse_heterodyne_filter_knee = self.get_config_option('heterodyne', 'filter_knee', 'float', default=0.25) self.coarse_heterodyne_sample_rate = self.get_config_option('heterodyne', 'coarse_sample_rate', 'int', default=16384) self.coarse_heterodyne_resample_rate = self.get_config_option('heterodyne', 'coarse_resample_rate', 'float', default=1.0) self.coarse_heterodyne_channels = self.get_config_option('heterodyne', 'channels', 'dict') self.coarse_heterodyne_binary_output = self.get_config_option('heterodyne', 'binary_output', 'boolean', default=True) self.coarse_heterodyne_gzip_output = self.get_config_option('heterodyne', 'gzip_coarse_output', 'boolean', default=False) self.coarse_heterodyne_request_memory = self.get_config_option('heterodyne', 'coarse_request_memory', 'int', allownone=True) self.coarse_heterodyne_max_data_length = self.get_config_option('heterodyne', 'coarse_max_data_length', 'int', default=512) if self.error_code != 0: return if self.coarse_heterodyne_binary_output and self.coarse_heterodyne_gzip_output: print("Warning... cannot output coarse heterdyned data as gzip and binary. Defaulting to binary output.") self.coarse_heterodyne_gzip_output = False self.fine_heterodyne_resample_rate = self.get_config_option('heterodyne', 'fine_resample_rate', 'string', default='1/60') self.fine_heterodyne_stddev_thresh = self.get_config_option('heterodyne', 'stddev_thresh', 'float', default=3.5) self.fine_heterodyne_gzip_output = self.get_config_option('heterodyne', 'gzip_fine_output', 'boolean', default=True) self.fine_heterodyne_request_memory = self.get_config_option('heterodyne', 'fine_request_memory', 'int', allownone=True) if self.error_code != 0: return # create heterodyne job chetjob = heterodyneJob(self.heterodyne_exec, univ=self.heterodyne_universe, accgroup=self.accounting_group, accuser=self.accounting_group_user, logdir=self.log_dir, rundir=self.run_dir, subprefix='coarse', requestmemory=self.coarse_heterodyne_request_memory) fhetjob = heterodyneJob(self.heterodyne_exec, univ=self.heterodyne_universe, accgroup=self.accounting_group, accuser=self.accounting_group_user, logdir=self.log_dir, rundir=self.run_dir, subprefix='fine', requestmemory=self.fine_heterodyne_request_memory) self.processed_files = {} # dictionary of processed (fine heterodyned files) self.fine_heterodyne_nodes = {} # dictionary for fine heterodyne nodes to use as parents to later jobs self.modified_pulsars_segment_list_tmp = {} getmodsciencesegs = {} # flag for modified pulsars setting whether script needs to get the science segment list segfiletimes = {} # (for unmodified pulsars) store segment files that have already been generated (keyed by the start time) to save repeated segment list finding calls for ifo in self.ifos: getmodsciencesegs[ifo] = True segfiletimes[ifo] = {} # check if channels is a single value for each IFO or a list for ifo in self.ifos: if ifo not in self.coarse_heterodyne_channels: print("Error... could channel not specified for '{}'".format(ifo), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return else: if isinstance(self.coarse_heterodyne_channels[ifo], string_types): # convert to list self.coarse_heterodyne_channels[ifo] = [self.coarse_heterodyne_channels[ifo]] elif not isinstance(self.coarse_heterodyne_channels[ifo], list): print("Error... channel must be a string or a list of strings", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return # loop through all pulsars for pname in self.analysed_pulsars: par = self.analysed_pulsars[pname] modified_pulsar = unmodified_pulsar = False # check if pulsar is new/has modified par file, or has an unmodified par file if par in self.modified_pulsars: modified_pulsar = True elif par in self.unmodified_pulsars: unmodified_pulsar = True # get ephemeris earthfile, sunfile, timefile = self.get_ephemeris(pppu.psr_par(par)) if self.error_code != 0: return segfiles = {} self.fine_heterodyne_nodes[pname] = {} self.processed_files[pname] = {} # loop through detectors and set up directory structure and get science segments on first pass for ifo in self.ifos[:]: # use self.ifos[:] to create a copy of the IFOs list, so that removal can be done if required https://stackoverflow.com/a/8544571/1862861 psrdir = os.path.join(self.preprocessing_base_dir[ifo], pname) self.mkdirs(psrdir) if self.error_code != 0: return datadir = os.path.join(psrdir, 'data') self.mkdirs(datadir) if self.error_code != 0: return coarsedir = os.path.join(datadir, 'coarse') self.mkdirs(coarsedir) if self.error_code != 0: return finedir = os.path.join(datadir, 'fine') self.mkdirs(finedir) if self.error_code != 0: return segfiles[ifo] = os.path.join(psrdir, 'segments.txt') # get segment list if modified_pulsar: if getmodsciencesegs[ifo]: self.modified_pulsars_segment_list_tmp[ifo] = os.path.join(self.preprocessing_base_dir[ifo], 'segments_tmp.txt') self.find_data_segments(self.initial_start[ifo], self.endtime[ifo], ifo, self.modified_pulsars_segment_list_tmp[ifo]) if self.error_code != 0: return # copy segment list into pulsar directories try: shutil.copy2(self.modified_pulsars_segment_list_tmp[ifo], segfiles[ifo]) except: print("Error... could not copy segment list into pulsar directory", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return if getmodsciencesegs[ifo]: getmodsciencesegs[ifo] = False # remove any "previous_segments.txt" file if it exists prevsegs = os.path.join(os.path.dirname(segfiles[ifo]), 'previous_segments.txt') if os.path.isfile(prevsegs): try: os.remove(prevsegs) except: print("Warning... previous segment list file '%s' could not be removed" % prevsegs) elif unmodified_pulsar: # append current segment file (if it exists) to file containing previous segments if os.path.isfile(segfiles[ifo]): # get filename for list containing all previous segments prevsegs = os.path.join(os.path.dirname(segfiles[ifo]), 'previous_segments.txt') if not os.path.isfile(prevsegs): shutil.copy2(segfiles[ifo], prevsegs) # copy current (unupdated) segfile into previous segments file if it does not already exist self.segment_file_update.append((segfiles[ifo], prevsegs)) # add to list of files to be concatenated together at the end (i.e. provided no errors have occurred) if self.autonomous: # if running in automatic mode make sure we don't miss data by using the end time in the previous segment # file as the new start time (rather than self.starttime) if os.path.isfile(segfiles[ifo]): # get end time from segment file p = sp.check_output("tail -1 " + segfiles[ifo], shell=True) if len(p) != 0: self.starttime[ifo] = [int(p.split()[1])] # add as a list for consistency else: print("Error... could not get end time out of previous segment file '%s'" % segfiles[ifo], file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return # check if a segment file for the given start time has already been created, and if so use that if self.starttime[ifo][0] in segfiletimes[ifo]: shutil.copy2(segfiletimes[ifo][self.starttime[ifo][0]], segfiles[ifo]) else: segfiletimes[ifo][self.starttime[ifo][0]] = segfiles[ifo] # create new segment file self.find_data_segments(self.starttime[ifo], self.endtime[ifo], ifo, segfiles[ifo]) if self.error_code != 0: return # check if generated data segment file contains any segments and not in autonomous mode if self.warning_code == KNOPE_WARNING_NO_SEGMENTS and not self.autonomous: self.ifos.remove(ifo) # remove IFO for which no segments exist if len(self.ifos) == 0: print("Error... no segments were available for any required IFOs", file=sys.stderr) self.error_code = KNOPE_ERROR_NO_SEGMENTS return else: self.warning_code = 0 continue self.fine_heterodyne_nodes[pname][ifo] = {} self.processed_files[pname][ifo] = {} # loop through frequency factors for analysis for freqfactor in self.freq_factors: if not freqfactor%1.: # for integers just output directory as e.g. 2f freqfacdircoarse = os.path.join(coarsedir, '%df' % int(freqfactor)) freqfacdirfine = os.path.join(finedir, '%df' % int(freqfactor)) else: # for non-integers use 2 d.p. for dir name freqfacdircoarse = os.path.join(coarsedir, '%.2ff' % freqfactor) freqfacdirfine = os.path.join(finedir, '%.2ff' % freqfactor) self.mkdirs(freqfacdircoarse) if self.error_code != 0: return self.mkdirs(freqfacdirfine) if self.error_code != 0: return # loop over number of "datasets" finetmpfiles = [] # file list to concatenate if there is more than one "dataset" if self.ndatasets[ifo] > 1: # create concatenation node for fine heterodyne fineoutput = os.path.join(freqfacdirfine, 'fine-%s-%d-%d.txt' % (ifo, int(self.starttime[ifo][0]), int(self.endtime[ifo][-1]))) if self.fine_heterodyne_gzip_output: fineoutput += '.gz' concatjob = concatJob(subfile=os.path.join(freqfacdirfine, 'concat.sub'), output=fineoutput, accgroup=self.accounting_group, accuser=self.accounting_group_user, logdir=self.log_dir) for k in range(self.ndatasets[ifo]): # create output files if self.warning_code != KNOPE_WARNING_NO_SEGMENTS: # only add coarse heterodyne if there is segment data if self.coarse_heterodyne_binary_output: coarseoutput = os.path.join(freqfacdircoarse, 'coarse-%s-%d-%d.bin' % (ifo, int(self.starttime[ifo][k]), int(self.endtime[ifo][k]))) else: coarseoutput = os.path.join(freqfacdircoarse, 'coarse-%s-%d-%d.txt' % (ifo, int(self.starttime[ifo][k]), int(self.endtime[ifo][k]))) # create coarse heterodyne node coarsenode = heterodyneNode(chetjob) # add data find parent to coarse job if self.datafind_nodes is not None: coarsenode.add_parent(self.datafind_nodes[ifo][k]) # set data, segment location, channel and output location coarsenode.set_data_file(self.cache_files[ifo][k]) coarsenode.set_max_data_length(self.coarse_heterodyne_max_data_length) coarsenode.set_seg_list(segfiles[ifo]) coarsenode.set_channel(self.coarse_heterodyne_channels[ifo][k]) coarsenode.set_output_file(coarseoutput) if self.coarse_heterodyne_binary_output: coarsenode.set_binoutput() if self.coarse_heterodyne_gzip_output: coarsenode.set_gzip_output() # set some coarse heterodyne info coarsenode.set_ifo(ifo) # detector coarsenode.set_het_flag(0) # perform coarse heterodyne coarsenode.set_pulsar(pname) # pulsar name coarsenode.set_param_file(par) # pulsar parameter file coarsenode.set_filter_knee(self.coarse_heterodyne_filter_knee) coarsenode.set_sample_rate(self.coarse_heterodyne_sample_rate) coarsenode.set_resample_rate(self.coarse_heterodyne_resample_rate) coarsenode.set_freq_factor(freqfactor) # multiplicative factor for the rotation frequency self.add_node(coarsenode) if self.warning_code != KNOPE_WARNING_NO_SEGMENTS: # only add fine heterodyne if there is segment data # create fine heterodyne node finenode = heterodyneNode(fhetjob) if self.ndatasets[ifo] > 1: # create concatnode if k == 0: concatnode = concatNode(concatjob) concatnode.add_parent(finenode) # add fine heterodyne as parent # add coarse parent to fine job finenode.add_parent(coarsenode) # create output files fineoutput = os.path.join(freqfacdirfine, 'fine-%s-%d-%d.txt' % (ifo, int(self.starttime[ifo][k]), int(self.endtime[ifo][k]))) finetmpfiles.append(fineoutput) # set data, segment location and output location if self.coarse_heterodyne_gzip_output: finenode.set_data_file(coarseoutput+'.gz') else: finenode.set_data_file(coarseoutput) if self.coarse_heterodyne_binary_output: finenode.set_bininput() finenode.set_seg_list(segfiles[ifo]) finenode.set_output_file(fineoutput) # set fine heterodyne info finenode.set_ifo(ifo) finenode.set_het_flag(1) # perform fine heterodyne finenode.set_pulsar(pname) # pulsar name finenode.set_param_file(par) # pulsar parameter file finenode.set_sample_rate(self.coarse_heterodyne_resample_rate) # use resample rate from coarse heterodyne finenode.set_resample_rate(self.fine_heterodyne_resample_rate) finenode.set_filter_knee(self.coarse_heterodyne_filter_knee) finenode.set_freq_factor(freqfactor) finenode.set_stddev_thresh(self.fine_heterodyne_stddev_thresh) finenode.set_ephem_earth_file(earthfile) finenode.set_ephem_sun_file(sunfile) finenode.set_ephem_time_file(timefile) if self.fine_heterodyne_gzip_output: finenode.set_gzip_output() finetmpfiles[-1] += '.gz' self.add_node(finenode) if self.ndatasets[ifo] > 1: if k == self.ndatasets[ifo]-1: concatnode.set_files(finetmpfiles) self.add_node(concatnode) # reset output name to that for the concatenated file fineoutput = os.path.join(freqfacdirfine, 'fine-%s-%d-%d.txt' % (ifo, int(self.starttime[ifo][0]), int(self.endtime[ifo][-1]))) # remove the seperate find heterodyned files rmjob = removeJob(accgroup=self.accounting_group, accuser=self.accounting_group_user, logdir=self.log_dir, rundir=self.run_dir) # job for removing old files rmnode = removeNode(rmjob) rmnode.set_files(finetmpfiles) rmnode.add_parent(concatnode) self.add_node(rmnode) # record fine nodes and files if self.ndatasets[ifo] > 1: self.fine_heterodyne_nodes[pname][ifo][freqfactor] = concatnode # set concatenation node as parent of future children else: self.fine_heterodyne_nodes[pname][ifo][freqfactor] = finenode if self.fine_heterodyne_gzip_output: fineoutput = fineoutput + '.gz' # append .gz to recorded file name if modified_pulsar: if self.warning_code != KNOPE_WARNING_NO_SEGMENTS: self.processed_files[pname][ifo][freqfactor] = [fineoutput] elif unmodified_pulsar: # check for other fine heterodyned files already in the directory hetfilescheck = [os.path.join(freqfacdirfine, hf) for hf in os.listdir(freqfacdirfine) if 'fine' in hf] hetfilescheck.sort() # if gzipping the current file, but older files are not gzipped, gzip them now if self.fine_heterodyne_gzip_output: for i, hf in enumerate(hetfilescheck): if '.gz' not in hf: # try gzipping file p = sp.Popen("gzip " + hf, shell=True) out, err = p.communicate() if p.returncode != 0: print("Error... could not gzip previous fine heterodyned file '%s': %s, %s" % (hf, out, err), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return hetfilescheck[i] = hf + '.gz' # add .gz suffix else: # alternatively, if not gzipping current file, but older files are gzipped then gunzip them now for i, hf in enumerate(hetfilescheck): if '.gz' in hf: # try gunzipping file p = sp.Popen("gunzip " + hf, shell=True) out, err = p.communicate() if p.returncode != 0: print("Error... could not gunzip previous fine heterodyned file '%s': %s, %s" % (hf, out, err), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return self.processed_files[pname][ifo][freqfactor] = hetfilescheck if self.warning_code == KNOPE_WARNING_NO_SEGMENTS: # if there are no segments this time, and also no previous fine heterodyned data then ignore this IFO if len(hetfilescheck) == 0: self.ifos.remove(ifo) # remove IFO for which no segments exist if len(self.ifos) == 0: print("Error... no segments were available for any required IFOs", file=sys.stderr) self.error_code = KNOPE_ERROR_NO_SEGMENTS return else: self.warning_code = 0 continue self.processed_files[pname][ifo][freqfactor].append(fineoutput) # append new file # remove temporary segment file if len(self.modified_pulsars) > 0: for ifo in self.ifos: try: if os.path.isfile(self.modified_pulsars_segment_list_tmp[ifo]): os.remove(self.modified_pulsars_segment_list_tmp[ifo]) except: print("Warning... could not remove temporary segment list '%s'" % self.modified_pulsars_segment_list_tmp[ifo]) def setup_splinter(self): """ Setup the Spectral Interpolation jobs/nodes. """ # NOTE: splinter runs on multiple pulsars at once, so requires a directory containing the necessary pulsars, # therefore for new/modified par files need to be linked to from one directory, and unmodifed par files need # to be linked to from another directory # Spectral interpolation output files have the format Splinter_PSRJNAME_DET e.g. Splinter_J0534+2200_H1 # get executable self.splinter_exec = self.get_config_option('splinter', 'splinter_exec', default='lalapps_SplInter') if self.error_code != 0: return self.splinter_modified_pars = {} self.splinter_unmodified_pars = {} # check file exists and is executable if not os.path.isfile(self.splinter_exec) or not os.access(self.splinter_exec, os.X_OK): print("Warning... 'splinter_exec' in '[splinter]' does not exist or is not an executable. Try finding code in path.") splexec = self.find_exec_file('lalapps_SplInter') if splexec == None: print("Error... could not find 'lalapps_SplInter' in 'PATH'", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return else: self.splinter_exec = splexec # Get Condor universe (default to vanilla) self.splinter_universe = self.get_config_option('splinter', 'universe', default='vanilla') if self.error_code != 0: return # check any memory requirements for a splinter job self.splinter_request_memory = self.get_config_option('splinter', 'splinter_request_memory', 'int', allownone=True) # create splinter job spljob = splinterJob(self.splinter_exec, univ=self.splinter_universe, accgroup=self.accounting_group, accuser=self.accounting_group_user, logdir=self.log_dir, rundir=self.run_dir, requestmemory=self.splinter_request_memory) # create job for moving splinter output files mvjob = moveJob(accgroup=self.accounting_group, accuser=self.accounting_group_user, logdir=self.log_dir, rundir=self.run_dir) # get splinter options self.splinter_bandwidth = self.get_config_option('splinter', 'bandwidth', 'float', 0.3) self.splinter_freq_range = self.get_config_option('splinter', 'freq_range', 'list', [30., 2000.]) self.splinter_stddev_thresh = self.get_config_option('splinter', 'stddev_thresh', 'float', 3.5) self.splinter_min_seg_length = self.get_config_option('splinter', 'min_seg_length', 'int', 1800) self.splinter_max_seg_length = self.get_config_option('splinter', 'max_seg_length', 'int', 21600) # max seg length defaults to 1/4 day self.splinter_gzip_output = self.get_config_option('splinter', 'gzip_output', 'boolean', False) if self.error_code != 0: return if self.splinter_min_seg_length > self.splinter_max_seg_length: print("Error... minimum segment length ({}) for SplInter is larger than maximum segment length ({})".format(self.splinter_min_seg_length, self.splinter_max_seg_length), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return # create splinter nodes (one for unmodified pulsars and one for modified pulsars) self.splinter_nodes_modified = {} self.splinter_nodes_unmodified = {} self.splinter_modified_pulsar_dir = None self.splinter_unmodified_pulsar_dir = None # set location for SplInter data to eventually be copied to self.splinter_data_location = {} self.processed_files = {} # dictionary of processed (splinter) files # check for whether there are any modified or unmodified pulsars modpars = False unmodpars = False if len(self.modified_pulsars) > 0: modpars = True if len(self.unmodified_pulsars) > 0: unmodpars = True for ifo in self.ifos[:]: # iterate over copy of self.ifos as may need to remove IFOs if no segments found # Create splinter base directory to contain data products and par file directories self.splinter_dir = os.path.join(self.preprocessing_base_dir[ifo], 'splinter') self.mkdirs(self.splinter_dir) if self.error_code != 0: return # Create par file directories: one containing links to files for modified pulsars and one for unmodified pulsars - if directories already exist then remove files from them for modsuffix in ['modified', 'unmodified']: pardir = os.path.join(self.splinter_dir, modsuffix) if modsuffix == 'modified': if modpars is True: self.splinter_modified_pulsar_dir = pardir parlist = self.modified_pulsars self.splinter_nodes_modified[ifo] = {} else: continue if modsuffix == 'unmodified': if unmodpars is True: self.splinter_unmodified_pulsar_dir = pardir parlist = self.unmodified_pulsars self.splinter_nodes_unmodified[ifo] = {} else: continue if os.path.isdir(pardir): # remove all files within the directory for f in os.listdir(pardir): try: os.remove(os.path.join(pardir, f)) except: print("Warning... problem removing par file '%s' from '%s'. This file may be overwritten." % (f, pardir)) else: # make the directory and make symbolic links to all the modified par files in to self.mkdirs(pardir) if self.error_code != 0: return for par in parlist: try: psr = pppu.psr_par(par) pname = psr['PSRJ'] parlink = os.path.join(pardir, os.path.basename(par)) if modsuffix == 'modified': self.splinter_modified_pars[pname] = parlink if modsuffix == 'unmodified': self.splinter_unmodified_pars[pname] = parlink os.symlink(par, parlink) # create dictionary to hold list of names of the output files that will be created by lalapps_SplInter if pname not in self.processed_files: self.processed_files[pname] = {} self.splinter_data_location[pname] = {} self.processed_files[pname][ifo] = {} self.splinter_data_location[pname][ifo] = {} # create directory for the output file to eventually moved into psrdir = os.path.join(self.preprocessing_base_dir[ifo], pname) self.mkdirs(psrdir) if self.error_code != 0: return datadir = os.path.join(psrdir, 'data') self.mkdirs(datadir) if self.error_code != 0: return splintercpydir = os.path.join(datadir, 'splinter') self.mkdirs(splintercpydir) if self.error_code != 0: return for freqfactor in self.freq_factors: if not freqfactor%1.: # for integers just output director as e.g. 2f ffdir = os.path.join(splintercpydir, '%df' % int(freqfactor)) else: # for non-integers use 2 d.p. for dir name ffdir = os.path.join(splintercpydir, '%.3ff' % int(freqfactor)) # the name of the file that will be output by lalapps_SplInter self.processed_files[pname][ifo][freqfactor] = [os.path.join(ffdir, 'SplInter_%s_%s' % (pname, ifo))] if self.splinter_gzip_output: self.processed_files[pname][ifo][freqfactor][0] += '.gz' # add gz suffix for gzipped files # the location to move that file to self.mkdirs(ffdir) if self.error_code != 0: return self.splinter_data_location[pname][ifo][freqfactor] = ffdir except: print("Warning... could not create link to par file '%s' in '%s'. This file may be overwritten." % (par, pardir)) # reset modpars and unmodpars based on whether any files are in self.splinter_(un)modified_pars if len(self.splinter_modified_pars) == 0: modpars = False if len(self.splinter_unmodified_pars) == 0: unmodpars = False # Create segment list for modified pulsars modsegfile = None if modpars: modsegfile = os.path.join(self.preprocessing_base_dir[ifo], 'segments_modified.txt') self.find_data_segments(self.initial_start[ifo], self.endtime[ifo], ifo, modsegfile) if self.error_code != 0: return if self.warning_code == KNOPE_WARNING_NO_SEGMENTS: self.ifos.remove(ifo) # remove current IFO from job if len(self.ifos) == 0: print("Error... no segments could were available for any required IFOs", file=sys.stderr) self.error_code = KNOPE_ERROR_NO_SEGMENTS return else: self.warning_code = 0 continue # Create segment list for unmodified pulsars unmodsegfile = None if unmodpars: unmodsegfile = os.path.join(self.preprocessing_base_dir[ifo], 'segments_unmodified.txt') if self.autonomous: # if running in automatic mode make sure we don't miss data by using the end time in the previous segment # file as the new start time (rather than self.starttime) if os.path.isfile(unmodsegfile): # get end time from segment list p = sp.check_output("tail -1 " + unmodsegfile, shell=True) if len(p) != 0: self.starttime[ifo] = int(p.split()[1]) else: print("Error... could not get end time out of previous segment file '%s'" % unmodsegfile, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return # create new segment file self.find_data_segments(self.starttime[ifo], self.endtime[ifo], ifo, unmodsegfile) if self.error_code != 0: return if self.warning_code == KNOPE_WARNING_NO_SEGMENTS: self.ifos.remove(ifo) # remove current IFO from job if len(self.ifos) == 0: print("Error... no segments were available for any required IFOs", file=sys.stderr) self.error_code = KNOPE_ERROR_NO_SEGMENTS return else: self.warning_code = 0 continue # append segments to file containing all previously analysed segments # get filename for list containing all previous segments prevsegs = os.path.join(self.preprocessing_base_dir[ifo], 'previous_segments.txt') if not os.path.isfile(prevsegs): shutil.copy2(unmodsegfile, prevsegs) # copy current (unupdated) segfile into previous segments file if it does not already exist self.segment_file_update.append((unmodsegfile, prevsegs)) # add to list of files to be concatenated together at the end (i.e. provided no errors have occurred) # loop through frequency factors for analysis for freqfactor in self.freq_factors: if not freqfactor%1.: # for integers just output director as e.g. 2f splinterdir = os.path.join(self.splinter_dir, '%df' % int(freqfactor)) else: # for non-integers use 2 d.p. for dir name splinterdir = os.path.join(self.splinter_dir, '%.2ff' % freqfactor) self.mkdirs(splinterdir) if self.error_code != 0: return # create nodes (first for modified pulsars and second for unmodified) splsegfiles = [modsegfile, unmodsegfile] splpardirs = [self.splinter_modified_pulsar_dir, self.splinter_unmodified_pulsar_dir] for idx, splbool in enumerate([modpars, unmodpars]): if splbool: splnode = splinterNode(spljob) splnode.add_parent(self.datafind_nodes[ifo][0]) splnode.set_sft_lalcache(self.cache_files[ifo]) # SFT cache file splnode.set_output_dir(splinterdir) # output directory splnode.set_param_dir(splpardirs[idx]) # pulsar parameter file directory splnode.set_seg_list(splsegfiles[idx]) # science segment list file splnode.set_freq_factor(freqfactor) # frequency scaling factor splnode.set_ifo(ifo) # detector splnode.set_bandwidth(self.splinter_bandwidth) # bandwidth of data to use splnode.set_min_seg_length(self.splinter_min_seg_length) # minimum length of usable science segments splnode.set_max_seg_length(self.splinter_max_seg_length) # maximum segment length (use to stop too many SFTs being read in at once and eating all the memory) splnode.set_start_freq(self.splinter_freq_range[0]) # low frequency cut-off to read from SFTs splnode.set_end_freq(self.splinter_freq_range[1]) # high frequency cut-off to read from SFTs splnode.set_stddev_thresh(self.splinter_stddev_thresh) # threshold for vetoing data splnode.set_ephem_dir(self.ephem_path) # path to solar system ephemeris files if self.splinter_gzip_output: splnode.set_gzip_output() # gzip the output files if idx == 0: self.splinter_nodes_modified[ifo][freqfactor] = splnode else: self.splinter_nodes_unmodified[ifo][freqfactor] = splnode self.add_node(splnode) # add node to DAG # move all created files into the required directory structure for par in parlist: psr = pppu.psr_par(par) pname = psr['PSRJ'] mvnode = moveNode(mvjob) splinterfile = os.path.join(splinterdir, 'SplInter_%s_%s' % (pname, ifo)) if self.splinter_gzip_output: splinterfile += '.gz' mvnode.set_source(splinterfile) mvnode.set_destination(self.processed_files[pname][ifo][freqfactor][0]) mvnode.add_parent(splnode) self.add_node(mvnode) def get_ephemeris(self, psr): """ Get the ephemeris file information based on the content of the psr file """ # get ephemeris value from par file ephem = psr['EPHEM'] if ephem is None: # if not present default to DE405 ephem = 'DE405' # get the time units (TDB) units = psr['UNITS'] if units is None: # default to TEMPO2 standard of TCB units = 'TCB' earthfile = os.path.join(self.ephem_path, 'earth00-40-%s.dat.gz' % ephem) sunfile = os.path.join(self.ephem_path, 'sun00-40-%s.dat.gz' % ephem) if units == 'TDB': timefile = os.path.join(self.ephem_path, 'tdb_2000-2040.dat.gz') else: timefile = os.path.join(self.ephem_path, 'te405_2000-2040.dat.gz') if not os.path.isfile(earthfile): print("Error... Earth ephemeris file '%s' does not exist" % earthfile, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return if not os.path.isfile(sunfile): print("Error... Sun ephemeris file '%s' does not exist" % sunfile, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return if not os.path.isfile(timefile): print("Error... time ephemeris file '%s' does not exist" % timefile, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return return earthfile, sunfile, timefile def concatenate_files(self): """ Create jobs to concatenate fine heterodyned/Splinter data files in cases where multiple files exist (e.g. in automated search). Go though the processed file list and find ones that have more than one file, also also remove previous files. """ rmjob = removeJob(accgroup=self.accounting_group, accuser=self.accounting_group_user, logdir=self.log_dir, rundir=self.run_dir) # job for removing old files self.remove_job = rmjob self.concat_nodes = {} # concatenation job nodes # loop over pulsars for pitem in self.processed_files: # pitems will be pulsar names pifos = self.processed_files[pitem] # loop through detectors for pifo in pifos: ffs = pifos[pifo] # loop through frequency factors for ff in ffs: finefiles = ffs[ff] concatnode = None subloc = None prevfile = None newfinefile = None catfiles = None # check for more than one file if self.engine == 'heterodyne': # do this for heterodyned data if len(finefiles) > 1: if pitem not in self.concat_nodes: self.concat_nodes[pitem] = {} if pifo not in self.concat_nodes[pitem]: self.concat_nodes[pitem][pifo] = {} # get start time of first file and end time of last file (assuming files have name fine-det-start-end.txt) firststart = os.path.basename(finefiles[0]).split('-')[2] lastend = os.path.basename(finefiles[-1]).split('-')[3] lastend = lastend.replace('.txt', '') lastend = lastend.replace('.gz', '') # create new file name for concatenated file newfinefile = os.path.join(os.path.dirname(finefiles[0]), 'fine-%s-%s-%s.txt' % (pifo, firststart, lastend)) if self.fine_heterodyne_gzip_output: newfinefile = newfinefile + '.gz' catfiles = finefiles try: # fine heterodyne node may not exist if (in autonomous mode) no new segments were found, but previous data did exist parent = self.fine_heterodyne_nodes[pitem][pifo][ff] # add equivalent fine heterodyned node as parent except: parent = None # create new concat job for each case (each one uses it's own sub file) subloc = os.path.join(os.path.dirname(finefiles[0]), 'concat.sub') elif self.engine == 'splinter': # do this for splinter data # check if final location already contains a file prevfile = os.listdir(self.splinter_data_location[pitem][pifo][ff]) if len(prevfile) > 1: print("Error... more than one previous Splinter file in directory '%s'" % self.splinter_data_location[pitem][pifo][ff]) self.error_code = KNOPE_ERROR_GENERAL return if pitem in self.splinter_modified_pars: # for modified pulsars set start time (for file name) from run values startname = str(self.starttime[pifo][0]) else: # a previous file exists use the start time from that (for file name) if len(prevfile) == 1: try: startname = os.path.basename(prevfile).split('-')[1] except: print("Warning... could not get previous start time from Splinter file name '%s'. Using start time from this run: '%d'" % (prevfile, self.starttime[pifo])) startname = str(self.starttime[pifo][0]) catfiles = [os.path.join(self.splinter_data_location[pitem][pifo][ff], prevfile[0]), finefiles] else: startname = str(self.starttime[pifo][0]) try: # splinter node may not exist if (in autonomous mode) no new segments were found, but previous data did exist parent = self.splinter_nodes_unmodified[pifo][ff] except: parent = None endname = str(self.endtime[pifo][-1]) # create new file name containing analysis time stamps newfinefile = os.path.join(self.splinter_data_location[pitem][pifo][ff], '%s-%s-%s.txt' % ((os.path.basename(finefiles[0])).strip('.gz'), startname, endname)) if self.splinter_gzip_output: newfinefile = newfinefile + '.gz' # create new concat job for each case (each one uses it's own sub file) subloc = os.path.join(self.splinter_data_location[pitem][pifo][ff], 'concat.sub') # create new concat job for each case (each one uses it's own sub file) if subloc is not None and catfiles is not None: concatjob = concatJob(subfile=subloc, output=newfinefile, accgroup=self.accounting_group, accuser=self.accounting_group_user, logdir=self.log_dir) concatnode = concatNode(concatjob) concatnode.set_files(catfiles) if parent is not None: concatnode.add_parent(parent) self.add_node(concatnode) self.concat_nodes[pitem][pifo][ff] = concatnode # make the processed file now just contain the final file location self.processed_files[pitem][pifo][ff] = [newfinefile] # remove original files if (len(finefiles) > 1 or self.engine == 'splinter') and concatnode != None: rmnode = removeNode(rmjob) if prevfile != None: # remove previous files finefiles.append(os.path.join(self.splinter_data_location[pitem][pifo][ff], prevfile[0])) rmnode.set_files(finefiles) if subloc is not None and catfiles is not None: rmnode.add_parent(concatnode) # only do after concatenation self.add_node(rmnode) def get_config_option(self, section, option, cftype=None, default=None, allownone=False): """ Get a value of type cftype ('string', 'int', 'float', 'boolean', 'list', 'dict' or 'dir') from the configuration parser object. Return value on success and None on failure """ value = None # return value for failure to parse option if cftype == None or cftype == 'string' or cftype == 'dir': try: value = self.config.get(section, option) except: if not allownone: if not isinstance(default, str): print("Error... could not parse '%s' option from '[%s]' section." % (option, section), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL else: print("Warning... could not parse '%s' option from '[%s]' section. Defaulting to %s." % (option, section, default)) value = default elif cftype == 'float': try: value = self.config.getfloat(section, option) except: if not allownone: if not isinstance(default, float): print("Error... could not parse '%s' float option from '[%s]' section." % (option, section), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL else: print("Warning... could not parse '%s' float option from '[%s]' section. Defaulting to %f." % (option, section, default)) value = default elif cftype == 'boolean': try: value = self.config.getboolean(section, option) except: if not allownone: if not isinstance(default, bool): print("Error... could not parse '%s' boolean option from '[%s]' section." % (option, section), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL else: print("Warning... could not parse '%s' boolean option from '[%s]' section. Defaulting to %r." % (option, section, default)) value = default elif cftype == 'int': try: value = self.config.getint(section, option) except: if not allownone: if not isinstance(default, int): print("Error... could not parse '%s' int option from '[%s]' section." % (option, section), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL else: print("Warning... could not parse '%s' int option from '[%s]' section. Defaulting to %d." % (option, section, default)) value = default elif cftype == 'list': try: value = ast.literal_eval(self.config.get(section, option)) if not isinstance(value, list): value = [value] except: if not allownone: if not isinstance(default, list): print("Error... could not parse '%s' list option from '[%s]' section." % (option, section), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL else: print("Warning... could not parse '%s' list option from '[%s]' section. Defaulting to [%s]." % (option, section, ', '.join(default))) value = default elif cftype == 'dict': try: value = ast.literal_eval(self.config.get(section, option)) if not isinstance(value, dict) and isinstance(default, dict): print("Warning... could not parse '%s' dictionary option from '[%s]' section. Defaulting to %s." % (option, section, str(default))) value = default elif not isinstance(value, dict) and not isinstance(default, dict): if not allownone: print("Error... could not parse '%s' dictionary option from '[%s]' section." % (option, section), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL else: value = None except: if not allownone: if not isinstance(default, dict): print("Error... could not parse '%s' dictionary option from '[%s]' section." % (option, section), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL else: print("Warning... could not parse '%s' dictionary option from '[%s]' section. Defaulting to %s." % (option, section, str(default))) value = default else: print("Error... unknown trying to get unknown type '%s' from configuration file" % cftype, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return value def check_universe(self, universe): """ Check Condor universe is 'local', 'vanilla' or 'standard' """ if universe not in ['local', 'vanilla', 'standard']: return True else: return False def setup_datafind(self): """ Create and setup the data find jobs. """ # Create data find job self.cache_files = {} self.datafind_job = None if self.config.has_option('condor', 'datafind'): # check if a dictionary of ready made cache files has been given datafind = self.get_config_option('condor', 'datafind', cftype='dict', allownone=True) if datafind is not None: # check there is a file for each detector and that they exist for ifo in self.ifos: if ifo not in datafind: print("Warning... no frame/SFT cache file given for %s, try using system gw_data_find instead" % ifo) datafindexec = self.find_exec_file('gw_data_find') if datafindexec is None: print("Error... could not find 'gw_data_find' in your 'PATH'", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return else: self.config.set('condor', 'datafind', datafindexec) # set value in config file parser else: if not isinstance(datafind[ifo], list): datafind[ifo] = [datafind[ifo]] self.cache_files[ifo] = [] for cachefile in datafind[ifo]: if not os.path.isfile(cachefile): print("Warning... frame/SFT cache file '%s' does not exist, try using system gw_data_find instead" % ifo) datafindexec = self.find_exec_file('gw_data_find') if datafindexec is None: print("Error... could not find 'gw_data_find' in your 'PATH'", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return else: self.config.set('condor', 'datafind', datafindexec) # set value in config file parser else: self.cache_files[ifo].append(cachefile) # check if a datafind job is needed for any of the detectors if len(self.cache_files) < len(self.ifos): self.datafind_job = pipeline.LSCDataFindJob(list(self.preprocessing_base_dir.values())[0], self.log_dir, self.config) else: # a data find exectable has been given datafind = self.get_config_option('condor', 'datafind') if os.path.isfile(datafind) and os.access(datafind, os.X_OK): self.datafind_job = pipeline.LSCDataFindJob(list(self.preprocessing_base_dir.values())[0], self.log_dir, self.config) else: print("Warning... data find executable '%s' does not exist, or is not executable, try using system gw_data_find instead" % datafind) datafindexec = self.find_exec_file('gw_data_find') if datafindexec is None: print("Error... could not find 'gw_data_find' in your 'PATH'", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return else: self.config.set('condor', 'datafind', datafindexec) # set value in config file parser self.datafind_job = pipeline.LSCDataFindJob(list(self.preprocessing_base_dir.values())[0], self.log_dir, self.config) else: # if no data find is specified try using the system gw_data_find datafindexec = self.find_exec_file('gw_data_find') if datafindexec is None: print("Error... could not find 'gw_data_find' in your 'PATH'", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return else: self.config.set('condor', 'datafind', datafindexec) # set value in config file parser self.datafind_job = pipeline.LSCDataFindJob(list(self.preprocessing_base_dir.values())[0], self.log_dir, self.config) # add additional options to data find job if self.datafind_job is not None: self.datafind_job.add_condor_cmd('accounting_group', self.accounting_group) if self.accounting_group_user is not None: self.datafind_job.add_condor_cmd('accounting_group_user', self.accounting_group_user) # reset the sub file location self.datafind_job.set_sub_file(os.path.join(self.run_dir, 'datafind.sub')) # Set gw_data_find nodes (one for each detector) if len(self.cache_files) < len(self.ifos): self.set_datafind_nodes() else: self.datafind_nodes = None def set_datafind_nodes(self): """ Add data find nodes to a dictionary of nodes """ self.datafind_nodes = {} if self.config.has_option('datafind', 'type'): frtypes = ast.literal_eval(self.config.get('datafind', 'type')) if not isinstance(frtypes, dict): print("Error... the data find 'types' must be a dictionary of values for each detector", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return else: print("Error... no frame 'type' specified for data find", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return for ifo in self.ifos: frtypelist = [] # list of frame types if not ifo in frtypes: print("Error... no data find type for %s" % ifo, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return else: if isinstance(frtypes[ifo], string_types): for i in range(self.ndatasets[ifo]): frtypelist.append(frtypes[ifo]) elif isinstance(frtypes[ifo], list): frtypelist = frtypes[ifo] if len(frtypelist) != self.ndatasets[ifo]: print("Error... the number of frame types must be the same as the number of start times", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return else: print("Error... frame types for '{}' must be a single string or a list of strings".format(ifo), file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return # check if cache file is already set for the given detector if ifo in self.cache_files: continue self.datafind_nodes[ifo] = [] self.cache_files[ifo] = [] if self.engine == 'splinter': # for splinter we can concatenate the data find output SFT caches together concatjob = concatJob(subfile=os.path.join(self.preprocessing_base_dir[ifo], 'concat.sub'), output=os.path.join(self.preprocessing_base_dir[ifo], 'cache.lcf'), accgroup=self.accounting_group, accuser=self.accounting_group_user, logdir=self.log_dir) for i in range(self.ndatasets[ifo]): dfnode = pipeline.LSCDataFindNode(self.datafind_job) dfnode.set_observatory(ifo[0]) dfnode.set_type(frtypelist[i]) dfnode.set_start(self.initial_start[ifo][i]) dfnode.set_end(self.endtime[ifo][i]) # reset the default LSCDataFindJob output cache filename if self.ndatasets[ifo] == 1: cachefile = os.path.join(self.preprocessing_base_dir[ifo], 'cache.lcf') else: cachefile = os.path.join(self.preprocessing_base_dir[ifo], 'cache_{}.lcf'.format(i)) self.cache_files[ifo].append(cachefile) dfnode.set_output(cachefile) self.datafind_nodes[ifo].append(dfnode) self.add_node(dfnode) # add nodes into DAG if self.engine == 'splinter': # for splinter we can concatenate the data find output SFT caches together if self.ndatasets[ifo] == 0: concatnode = concatNode(concatjob) concatnode.set_files(self.cache_files[ifo]) for dfn in self.datafind_nodes[ifo]: concatnode.add_parent(dfn) self.add_node(concatnode) # reset the self.datafind_nodes to the concatnode for future parents to use self.datafind_nodes[ifo] = [concatnode] # reset name of cache file self.cache_files[ifo] = [os.path.join(self.preprocessing_base_dir[ifo], 'cache.lcf')] def find_data_segments(self, starttime, endtime, ifo, outfile): """ Find data segments and output them to an acsii text segment file containing the start and end times of each segment. Pairs of start and end times in the `starttime` and `endtime` lists will be iterated over and all segments within those pair will be concatenated into the final file. Args: starttime (list): a list of GPS start times endtime (list): a list of GPS end times ifo (str): a detector name, e.g., 'H1' outfile (str): the file to output the segment list file to """ # check if segment file(s) given (as a dictionary) segfiles = self.get_config_option('segmentfind', 'segfind', cftype='dict', allownone=True) if segfiles is not None: # check if it is a dictionary if ifo not in segfiles: print("Error... No segment file given for '%s'" % ifo) self.error_code = KNOPE_ERROR_GENERAL return else: segfile = segfiles[ifo] # check segment file exists if not os.path.isfile(segfile): print("Error... segment file '%s' does not exist." % segfile) self.error_code = KNOPE_ERROR_GENERAL return else: # copy segment file to the 'outfile' location try: shutil.copyfile(segfile, outfile) except: print("Error... could not copy segment file to location of '%s'." % outfile) self.error_code = KNOPE_ERROR_GENERAL return # exit function # otherwise try and get the segment list using gwpy try: from gwpy.segments import DataQualityFlag except ImportError: print("Error... gwpy is required for finding segment lists") self.error_code = KNOPE_ERROR_GENERAL return # get server server = self.get_config_option('segmentfind', 'server', default='https://segments.ligo.org') # get segment types to include segmenttypes = self.get_config_option('segmentfind', 'segmenttype', cftype='dict') # get segment types to exclude excludesegs = self.get_config_option('segmentfind', 'excludetype', cftype='dict', allownone=True) # check segment types dictionary contains type for the given detector if ifo not in segmenttypes: print("Error... No segment type for %s" % ifo, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return if isinstance(starttime, int) and isinstance(endtime, int): sts = [starttime] # have start time as a list ets = [endtime] # have end time as a list elif isinstance(starttime, list) and isinstance(endtime, list): if len(starttime) != len(endtime): print("Error... list of start and end times for the segments are not the same lengths", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return for st, et in zip(starttime, endtime): if st > et: print("Error... start time comes before end time!", file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL return sts = starttime ets = endtime else: print("Error... start and end times must be integers or lists of integers") self.error_code = KNOPE_ERROR_GENERAL return sidx = 0 try: segfp = open(outfile, 'w') except IOError: print("Error... could not open segment list file") self.error_code = KNOPE_ERROR_GENERAL return # loop over start and end time pairs for st, et in zip(sts, ets): # check whether there are different segment types for the different times if isinstance(segmenttypes[ifo], list): if not isinstance(segmenttypes[ifo][sidx], string_types): print("Error... segment types must be a string") self.error_code = KNOPE_ERROR_GENERAL return elif len(segmenttypes[ifo]) != len(sts): print("Error... number of segment types is not the same as the number of start times") self.error_code = KNOPE_ERROR_GENERAL return else: segmenttype = segmenttypes[ifo][sidx] else: if not isinstance(segmenttypes[ifo], string_types): print("Error... segment types must be a string") self.error_code = KNOPE_ERROR_GENERAL return else: segmenttype = segmenttypes[ifo] # is multiple segment types are specified (comma seperated string), loop # over them and find the intersection segtypes = [sts.strip() for sts in segmenttype.split(',')] # split and strip whitespace segquery = None for i in range(len(segtypes)): # create query query = DataQualityFlag.query_dqsegdb(segtypes[i], st, et, url=server) query = query.active # use "active" segment list if excludesegs is not None: # check whether there are different exclusion types for the different times if ifo in excludesegs: if isinstance(excludesegs[ifo], list): if not isinstance(excludesegs[ifo][sidx], string_types): print("Error... exclude types must be a string") self.error_code = KNOPE_ERROR_GENERAL return elif len(excludesegs[ifo]) != len(sts): print("Error... number of exclude types is not the same as the number of start times") self.error_code = KNOPE_ERROR_GENERAL return else: excludetype = excludesegs[ifo][sidx] else: if not isinstance(excludesegs[ifo], string_types): print("Error... exclude types must be a string") self.error_code = KNOPE_ERROR_GENERAL return else: excludetype = excludesegs[ifo] if len(excludetype) > 0: segextypes = [sts.strip() for sts in excludetype.split(',')] # split and strip whitespace for j in range(len(segextypes)): exquery = DataQualityFlag.query_dqsegdb(segextypes[j], st, et, url=server) # exclude segments query = query & ~exquery.active if segquery is None: segquery = query.copy() else: # get intersection of segments segquery = segquery & query # write out segment list for thisseg in segquery: print(int(thisseg[0]), int(thisseg[1]), file=segfp) sidx += 1 segfp.close() def find_exec_file(self, filename): """ Search through the PATH environment for the first instance of the executable file "filename" """ if os.path.isfile(filename) and os.access(filename, os.X_OK): return filename # already have the executable file # search through PATH environment variable for d in os.environ['PATH'].split(os.pathsep): filecheck = os.path.join(d, filename) if os.path.isfile(filecheck) and os.access(filecheck, os.X_OK): return filecheck return None def mkdirs(self, path): """ Helper function. Make the given directory, creating intermediate dirs if necessary, and don't complain about it already existing. """ if os.access(path,os.W_OK) and os.path.isdir(path): return else: try: os.makedirs(path) except: print("Error... cannot make directory '%s'" % path, file=sys.stderr) self.error_code = KNOPE_ERROR_GENERAL class heterodyneJob(pipeline.CondorDAGJob, pipeline.AnalysisJob): """ A lalapps_heterodyne_pulsar job to heterodyne the data. """ def __init__(self, execu, univ='vanilla', accgroup=None, accuser=None, logdir=None, rundir=None, subprefix='', requestmemory=None): self.__executable = execu self.__universe = univ pipeline.CondorDAGJob.__init__(self, self.__universe, self.__executable) pipeline.AnalysisJob.__init__(self, None) if accgroup != None: self.add_condor_cmd('accounting_group', accgroup) if accuser != None: self.add_condor_cmd('accounting_group_user', accuser) self.add_condor_cmd('getenv','True') if requestmemory is not None: if isinstance(requestmemory, int): self.add_condor_cmd('request_memory', requestmemory) # set log files for job if logdir != None: self.set_stdout_file(os.path.join(logdir, 'heterodyne-$(cluster).out')) self.set_stderr_file(os.path.join(logdir, 'heterodyne-$(cluster).err')) else: self.set_stdout_file('heterodyne-$(cluster).out') self.set_stderr_file('heterodyne-$(cluster).err') if rundir != None: self.set_sub_file(os.path.join(rundir, subprefix+'heterodyne.sub')) else: self.set_sub_file(subprefix+'heterodyne.sub') class heterodyneNode(pipeline.CondorDAGNode, pipeline.AnalysisNode): """ A heterodyneNode runs an instance of lalapps_heterodyne_pulsar in a condor DAG. """ def __init__(self,job): """ job = A CondorDAGJob that can run an instance of lalapps_heterodyne_pulsar """ pipeline.CondorDAGNode.__init__(self,job) pipeline.AnalysisNode.__init__(self) # initilise job variables self.__ifo = None self.__param_file = None self.__freq_factor = None self.__filter_knee = None self.__sample_rate = None self.__resample_rate = None self.__data_file = None self.__channel = None self.__seg_list = None self.__data_file = None self.__output_file = None self.__het_flag = None self.__pulsar = None self.__high_pass = None self.__scale_fac = None self.__manual_epoch = None self.__gzip = False self.__max_data_length = None def set_data_file(self,data_file): # set file containing data to be heterodyne (either list of frames or coarse het output) self.add_var_opt('data-file',data_file) self.__data_file = data_file def set_max_data_length(self, maxdatalength): # set the maximum length of data to be read from a frame file self.add_var_opt('data-chunk-length', maxdatalength) self.__max_data_length = maxdatalength def set_output_file(self,output_file): # set output directory self.add_var_opt('output-file',output_file) self.__output_file = output_file def set_seg_list(self,seg_list): # set segment list to be used self.add_var_opt('seg-file',seg_list) self.__seg_list = seg_list def set_ifo(self, ifo): # set detector self.add_var_opt('ifo', ifo) self.__ifo = ifo def set_param_file(self, param_file): # set pulsar parameter file self.add_var_opt('param-file', param_file) self.__param_file = param_file def set_freq_factor(self, freq_factor): # set the factor by which to muliply the pulsar spin frequency (normally 2.0) self.add_var_opt('freq-factor', freq_factor) self.__freq_factor = freq_factor def set_param_file_update(self, param_file_update): # set file containing updated pulsar parameters self.add_var_opt('param-file-update',param_file_update) def set_manual_epoch(self, manual_epoch): # set manual pulsar epoch self.add_var_opt('manual-epoch',manual_epoch) self.__manual_epoch = manual_epoch def set_ephem_earth_file(self, ephem_earth_file): # set the file containing the earth's ephemeris self.add_var_opt('ephem-earth-file', ephem_earth_file) def set_ephem_sun_file(self, ephem_sun_file): # set the file containing the sun's ephemeris self.add_var_opt('ephem-sun-file', ephem_sun_file) def set_ephem_time_file(self, ephem_time_file): # set the file containing the Einstein delay correction ephemeris self.add_var_opt('ephem-time-file', ephem_time_file) def set_pulsar(self,pulsar): # set pulsar name self.add_var_opt('pulsar',pulsar) self.__pulsar = pulsar def set_het_flag(self,het_flag): # set heterodyne flag self.add_var_opt('heterodyne-flag',het_flag) self.__het_flag = het_flag def set_filter_knee(self,filter_knee): # set filter knee frequency self.add_var_opt('filter-knee',filter_knee) self.__filter_knee = filter_knee def set_channel(self,channel): # set channel containing data from frames self.add_var_opt('channel',channel) self.__channel = channel def set_sample_rate(self,sample_rate): # set sample rate of input data self.add_var_opt('sample-rate',sample_rate) self.__sample_rate = sample_rate def set_resample_rate(self,resample_rate): # set resample rate for output data self.add_var_opt('resample-rate',resample_rate) self.__resample_rate = resample_rate def set_stddev_thresh(self,stddev_thresh): # set standard deviation threshold at which to remove outliers self.add_var_opt('stddev-thresh',stddev_thresh) def set_calibrate(self): # set calibration flag self.add_var_opt('calibrate', '') # no variable required def set_verbose(self): # set verbose flag self.add_var_opt('verbose', '') # no variable required def set_bininput(self): # set binary input file flag self.add_var_opt('binary-input', '') # no variable required def set_binoutput(self): # set binary output file flag self.add_var_opt('binary-output', '') # no variable required def set_gzip_output(self): # set gzipped output file flag self.add_var_opt('gzip-output', '') # no variable required self.__gzip = True def set_response_function(self,response_function): # set reponse function file self.add_var_opt('response-file',response_function) def set_coefficient_file(self,coefficient_file): # set the file containing the calibration coefficients (e.g alpha and gammas) self.add_var_opt('coefficient-file',coefficient_file) def set_sensing_function(self,sensing_function): # set file containing the sensing function for calibration self.add_var_opt('sensing-function',sensing_function) def set_open_loop_gain(self,open_loop_gain): # set file containing the open loop gain for calibration self.add_var_opt('open-loop-gain',open_loop_gain) def set_scale_fac(self,scale_fac): # set scale factor for calibrated data self.add_var_opt('scale-factor',scale_fac) def set_high_pass(self,high_pass): # set high-pass frequency for calibrated data self.add_var_opt('high-pass-freq',high_pass) class splinterJob(pipeline.CondorDAGJob, pipeline.AnalysisJob): """ A lalapps_SplInter job to process SFT data. """ def __init__(self, execu, univ='vanilla', accgroup=None, accuser=None, logdir=None, rundir=None, requestmemory=None): self.__executable = execu self.__universe = univ pipeline.CondorDAGJob.__init__(self, self.__universe, self.__executable) pipeline.AnalysisJob.__init__(self, None) if accgroup != None: self.add_condor_cmd('accounting_group', accgroup) if accuser != None: self.add_condor_cmd('accounting_group_user', accuser) self.add_condor_cmd('getenv','True') if requestmemory is not None: if isinstance(requestmemory, int): self.add_condor_cmd('request_memory', requestmemory) # set log files for job if logdir != None: self.set_stdout_file(os.path.join(logdir, 'splinter-$(cluster).out')) self.set_stderr_file(os.path.join(logdir, 'splinter-$(cluster).err')) else: self.set_stdout_file('splinter-$(cluster).out') self.set_stderr_file('splinter-$(cluster).err') if rundir != None: self.set_sub_file(os.path.join(rundir, 'splinter.sub')) else: self.set_sub_file('splinter.sub') class splinterNode(pipeline.CondorDAGNode, pipeline.AnalysisNode): """ A splinterNode runs an instance of lalapps_Splinter in a condor DAG. """ def __init__(self,job): """ job = A CondorDAGJob that can run an instance of lalapps_SplInter """ pipeline.CondorDAGNode.__init__(self,job) pipeline.AnalysisNode.__init__(self) # initilise job variables self.__ifo = None self.__param_file = None self.__param_dir = None self.__start_freq = None self.__end_freq = None self.__freq_factor = None self.__data_file = None self.__seg_list = None self.__sft_cache = None self.__sft_loc = None self.__sft_lalcache = None self.__output_dir = None self.__stddev_thresh = None self.__bandwidth = None self.__min_seg_length = None self.__max_seg_length = None self.__starttime = None self.__endtime = None self.__ephem_dir = None self.__gzip = False def set_start_freq(self, f): # set the start frequency of the SFTs self.add_var_opt('start-freq', f) self.__start_freq = f def set_end_freq(self, f): # set the end frequency of the SFTs self.add_var_opt('end-freq', f) def set_sft_cache(self, f): # set file containing list of SFTs self.add_var_opt('sft-cache', f) self.__sft_cache = f def set_sft_lalcache(self, f): # set file containing list of SFTs in LALCache format self.add_var_opt('sft-lalcache', f) self.__sft_lalcache = f def set_sft_loc(self, f): # set directory of files containing list of SFTs self.add_var_opt('sft-loc', f) self.__sft_loc = f def set_output_dir(self, f): # set output directory self.add_var_opt('output-dir', f) self.__output_dir = f def set_seg_list(self, f): # set segment list to be used self.add_var_opt('seg-file', f) self.__seg_list = f def set_ifo(self, ifo): # set detector self.add_var_opt('ifo', ifo) self.__ifo = ifo def set_param_file(self, f): # set pulsar parameter file self.add_var_opt('param-file', f) self.__param_file = f def set_param_dir(self, f): # set pulsar parameter directory self.add_var_opt('param-dir', f) self.__param_dir = f def set_freq_factor(self, f): # set the factor by which to muliply the pulsar spin frequency (normally 2.0) self.add_var_opt('freq-factor', f) self.__freq_factor = f def set_bandwidth(self, f): # set the bandwidth to use in the interpolation self.add_var_opt('bandwidth', f) self.__bandwidth = f def set_min_seg_length(self, f): # set the minimum science segment length to use self.add_var_opt('min-seg-length', f) self.__min_seg_length = f def set_max_seg_length(self, f): # set the maximum length of a science segment to use self.add_var_opt('max-seg-length', f) self.__max_seg_length = f def set_ephem_dir(self, f): # set the directory containing the solar system ephemeris files self.add_var_opt('ephem-dir', f) self.__ephem_dir = f def set_stddev_thresh(self, f): # set standard deviation threshold at which to remove outliers self.add_var_opt('stddevthresh', f) self.__stddev_thresh = f def set_gzip_output(self): # set gzipped output file flag self.add_var_opt('gzip', '') # no variable required self.__gzip = True def set_starttime(self, f): # set the start time of data to use self.add_var_opt('starttime', f) self.__starttime = f def set_endtime(self, f): # set the end time of data to use self.add_var_opt('endtime', f) self.__endtime = f """ Job for concatenating processed (heterodyned or spectrally interpolated) files """ class concatJob(pipeline.CondorDAGJob, pipeline.AnalysisJob): """ A concatenation job (using "cat" and output to stdout - a new job is needed for each pulsar) """ def __init__(self, subfile=None, output=None, accgroup=None, accuser=None, logdir=None): self.__executable = '/bin/cat' # use cat self.__universe = 'local' # use local as "cat" should not be compute intensive pipeline.CondorDAGJob.__init__(self, self.__universe, self.__executable) pipeline.AnalysisJob.__init__(self, None) if accgroup != None: self.add_condor_cmd('accounting_group', accgroup) if accuser != None: self.add_condor_cmd('accounting_group_user', accuser) if subfile == None: print("Error... Condor sub file required for concatenation job", file=sys.stderr) sys.exit(1) else: self.set_sub_file(subfile) if output == None: print("Error... output file required for concatenation job", file=sys.stderr) sys.exit(1) else: self.set_stdout_file(output) if logdir != None: self.set_stderr_file(os.path.join(logdir, 'concat-$(cluster).err')) else: self.set_stderr_file('concat-$(cluster).err') self.add_arg('$(macrocatfiles)') # macro for input files to be concatenated class concatNode(pipeline.CondorDAGNode, pipeline.AnalysisNode): """ A node for a concatJob """ def __init__(self,job): """ job = A CondorDAGJob that can run an instance of parameter estimation code. """ pipeline.CondorDAGNode.__init__(self,job) pipeline.AnalysisNode.__init__(self) self.__files = None def set_files(self, files): # set the files to be concatenated ("files" is a list of files to be concatenated) self.add_macro('macrocatfiles', " ".join(files)) self.__files = files """ Job for removing files """ class removeJob(pipeline.CondorDAGJob, pipeline.AnalysisJob): """ A remove job (using "rm" to remove files) """ def __init__(self, accgroup=None, accuser=None, logdir=None, rundir=None): self.__executable = "/bin/rm" # use "rm" self.__universe = "local" # rm should not be compute intensive, so use local universe pipeline.CondorDAGJob.__init__(self, self.__universe, self.__executable) pipeline.AnalysisJob.__init__(self, None) if accgroup != None: self.add_condor_cmd('accounting_group', accgroup) if accuser != None: self.add_condor_cmd('accounting_group_user', accuser) # set log files for job if logdir != None: self.set_stdout_file(os.path.join(logdir, 'rm-$(cluster).out')) self.set_stderr_file(os.path.join(logdir, 'rm-$(cluster).err')) else: self.set_stdout_file('rm-$(cluster).out') self.set_stderr_file('rm-$(cluster).err') self.add_arg('-f $(macrormfiles)') # macro for input files to be removed (use "-f" to force removal) if rundir != None: self.set_sub_file(os.path.join(rundir, 'rm.sub')) else: self.set_sub_file('rm.sub') class removeNode(pipeline.CondorDAGNode, pipeline.AnalysisNode): """ An instance of a removeJob in a condor DAG. """ def __init__(self,job): """ job = A CondorDAGJob that can run an instance of rm. """ pipeline.CondorDAGNode.__init__(self,job) pipeline.AnalysisNode.__init__(self) self.__files = None self.set_retry(1) # retry the node once def set_files(self, files): # set file(s) to be removed, where "files" is a list containing all files self.add_macro('macrormfiles', " ".join(files)) self.__files = files """ Job for moving files """ class moveJob(pipeline.CondorDAGJob, pipeline.AnalysisJob): """ A move job (using "mv" to move files) """ def __init__(self, accgroup=None, accuser=None, logdir=None, rundir=None): self.__executable = "/bin/mv" # use "mv" self.__universe = "local" # mv should not be compute intensive, so use local universe pipeline.CondorDAGJob.__init__(self, self.__universe, self.__executable) pipeline.AnalysisJob.__init__(self, None) if accgroup != None: self.add_condor_cmd('accounting_group', accgroup) if accuser != None: self.add_condor_cmd('accounting_group_user', accuser) # set log files for job if logdir != None: self.set_stdout_file(os.path.join(logdir, 'mv-$(cluster).out')) self.set_stderr_file(os.path.join(logdir, 'mv-$(cluster).err')) else: self.set_stdout_file('mv-$(cluster).out') self.set_stderr_file('mv-$(cluster).err') self.add_arg('$(macrosource) $(macrodestination)') # macro for source and destination files if rundir != None: self.set_sub_file(os.path.join(rundir, 'mv.sub')) else: self.set_sub_file('mv.sub') class moveNode(pipeline.CondorDAGNode, pipeline.AnalysisNode): """ An instance of a moveJob in a condor DAG. """ def __init__(self,job): """ job = A CondorDAGJob that can run an instance of mv. """ pipeline.CondorDAGNode.__init__(self,job) pipeline.AnalysisNode.__init__(self) self.__sourcefile = None self.__destinationfile = None self.set_retry(1) # retry the node once def set_source(self, sfile): # set file to be moved self.add_macro('macrosource', sfile) self.__sourcefiles = sfile def set_destination(self, dfile): # set destination of file to be moved self.add_macro('macrodestination', dfile) self.__destinationfile = dfile """ Job for copying files """ class copyJob(pipeline.CondorDAGJob, pipeline.AnalysisJob): """ A copy job (using "cp" to copy files) """ def __init__(self, accgroup=None, accuser=None, logdir=None, rundir=None): self.__executable = "/bin/cp" # use "cp" self.__universe = "local" # cp should not be compute intensive, so use local universe pipeline.CondorDAGJob.__init__(self, self.__universe, self.__executable) pipeline.AnalysisJob.__init__(self, None) if accgroup != None: self.add_condor_cmd('accounting_group', accgroup) if accuser != None: self.add_condor_cmd('accounting_group_user', accuser) # set log files for job if logdir != None: self.set_stdout_file(os.path.join(logdir, 'cp-$(cluster).out')) self.set_stderr_file(os.path.join(logdir, 'cp-$(cluster).err')) else: self.set_stdout_file('cp-$(cluster).out') self.set_stderr_file('cp-$(cluster).err') self.add_arg('$(macrosource) $(macrodestination)') # macro for source and destination files if rundir != None: self.set_sub_file(os.path.join(rundir, 'cp.sub')) else: self.set_sub_file('cp.sub') class copyNode(pipeline.CondorDAGNode, pipeline.AnalysisNode): """ An instance of a copyJob in a condor DAG. """ def __init__(self,job): """ job = A CondorDAGJob that can run an instance of mv. """ pipeline.CondorDAGNode.__init__(self,job) pipeline.AnalysisNode.__init__(self) self.__sourcefile = None self.__destinationfile = None self.set_retry(1) # retry the node once def set_source(self, sfile): # set file to be moved self.add_macro('macrosource', sfile) self.__sourcefiles = sfile def set_destination(self, dfile): # set destination of file to be moved self.add_macro('macrodestination', dfile) self.__destinationfile = dfile """ Pulsar parameter estimation pipeline utilities """ class ppeJob(pipeline.CondorDAGJob, pipeline.AnalysisJob): """ A parameter estimation job """ def __init__(self, execu, univ='vanilla', accgroup=None, accuser=None, logdir=None, rundir=None, requestmemory=None): self.__executable = execu self.__universe = univ pipeline.CondorDAGJob.__init__(self, self.__universe, self.__executable) pipeline.AnalysisJob.__init__(self, None) if accgroup != None: self.add_condor_cmd('accounting_group', accgroup) if accuser != None: self.add_condor_cmd('accounting_group_user', accuser) self.add_condor_cmd('getenv','True') if requestmemory is not None: if isinstance(requestmemory, int): self.add_condor_cmd('request_memory', requestmemory) # set log files for job if logdir != None: self.set_stdout_file(os.path.join(logdir, 'ppe$(logname)-$(cluster).out')) self.set_stderr_file(os.path.join(logdir, 'ppe$(logname)-$(cluster).err')) else: self.set_stdout_file('ppe$(logname)-$(cluster).out') self.set_stderr_file('ppe$(logname)-$(cluster).err') if rundir != None: self.set_sub_file(os.path.join(rundir, 'ppe.sub')) else: self.set_sub_file('ppe.sub') self.add_arg('$(macroargs)') # macro for additional command line arguments that will not alway be used class ppeNode(pipeline.CondorDAGNode, pipeline.AnalysisNode): """ A pes runs an instance of the parameter estimation code in a condor DAG. """ def __init__(self,job,psrname=''): """ job = A CondorDAGJob that can run an instance of parameter estimation code. """ pipeline.CondorDAGNode.__init__(self,job) pipeline.AnalysisNode.__init__(self) # initilise job variables self.__detectors = None self.__par_file = None self.__cor_file = None self.__input_files = None self.__outfile = None self.__chunk_min = None self.__chunk_max = None self.__psi_bins = None self.__time_bins = None self.__prior_file = None self.__ephem_earth = None self.__ephem_sun = None self.__ephem_time = None self.__harmonics = None self.__biaxial = False self.__gaussian_like = False self.__randomise = None self.__starttime = None self.__endtime = None self.__truncate_time = None self.__truncate_samples = None self.__truncate_fraction = None self.__veto_threshold = None self.__Nlive = None self.__Nmcmc = None self.__Nmcmcinitial = None self.__Nruns = None self.__tolerance = None self.__randomseed = None self.__use_roq = False self.__roq_ntraining = None self.__roq_tolerance = None self.__roq_uniform = False self.__roq_output_weights = None self.__roq_input_weights = None self.__temperature = None self.__ensmeble_walk = None self.__ensemble_stretch = None self.__diffev = None self.__inject_file = None self.__inject_output = None self.__fake_data = None self.__fake_psd = None self.__fake_starts = None self.__fake_lengths = None self.__fake_dt = None self.__scale_snr = None self.__sample_files = None self.__sample_nlives = None self.__prior_cell = None # legacy inputs self.__oldChunks = None self.__sourceModel = False self.__verbose = False # set macroargs to be empty self.add_macro('macroargs', '') # set unique log name (set to pulsar name) if len(psrname) > 0: self.add_macro('logname', '-'+psrname) else: self.add_macro('logname', '') def set_detectors(self,detectors): # set detectors self.add_var_opt('detectors',detectors) self.__detectors = detectors # at detector names to logname if 'logname' in self.get_opts(): curmacroval = self.get_opts()['logname'] else: curmacroval = '' curmacroval = curmacroval + '-' + detectors.replace(',', '') self.add_macro('logname', curmacroval) def set_verbose(self): # set to run code in verbose mode self.add_var_opt('verbose', '') self.__verbose = True def set_par_file(self,parfile): # set pulsar parameter file self.add_var_opt('par-file',parfile) self.__par_file = parfile def set_cor_file(self,corfile): # set pulsar parameter correlation matrix file # add this into the generic 'macroargs' macro as it is not a value that is always required if 'macroargs' in self.get_opts(): curmacroval = self.get_opts()['macroargs'] else: curmacroval = '' curmacroval = curmacroval + ' --cor-file ' + corfile self.add_macro('macroargs', curmacroval) # self.add_var_opt('cor-file',corfile) self.__cor_file = corfile def set_input_files(self, inputfiles): # set data files for analysing self.add_var_opt('input-files', inputfiles) self.__input_files = inputfiles def set_outfile(self, of): # set the output file self.add_var_opt('outfile', of) self.__outfile = of def set_chunk_min(self, cmin): # set the minimum chunk length self.add_var_opt('chunk-min',cmin) self.__chunk_min = cmin def set_chunk_max(self, cmax): # set the maximum chunk length self.add_var_opt('chunk-max', cmax) self.__chunk_max = cmax def set_start_time(self, starttime): # set the start time of data to use self.add_var_opt('start-time', starttime) self.__starttime = starttime def set_end_time(self, endtime): # set the start time of data to use self.add_var_opt('end-time', endtime) self.__endtime = endtime self.__truncate_time = endtime def set_truncate_time(self, trunc): # set the truncation time to use (same as end time) if self.__endtime is None: self.add_var_opt('truncate-time', trunc) self.__truncate_time = trunc self.__endtime = trunc def set_truncate_samples(self, trunc): # set the truncation based on number of samples self.add_var_opt('truncate-samples', trunc) self.__truncate_samples = trunc def set_truncate_fraction(self, trunc): # set truncation based on fraction of the data self.add_var_opt('truncate-fraction', trunc) self.__truncate_fraction = trunc def set_veto_threshold(self, veto): # set value above which to veto data samples self.add_var_opt('veto-threshold', veto) self.__veto_threshold = veto def set_psi_bins(self, pb): # set the number of bins in psi for the psi vs time lookup table self.add_var_opt('psi-bins', pb) self.__psi_bins = pb def set_time_bins(self, tb): # set the number of bins in time for the psi vs time lookup table self.add_var_opt('time-bins',tb) self.__time_bins = tb def set_prior_file(self,pf): # set the prior ranges file self.add_var_opt('prior-file',pf) self.__prior_file = pf def set_ephem_earth(self,ee): # set earth ephemeris file self.add_var_opt('ephem-earth',ee) self.__ephem_earth = ee def set_ephem_sun(self,es): # set sun ephemeris file self.add_var_opt('ephem-sun',es) self.__ephem_sun = es def set_ephem_time(self,et): # set time correction ephemeris file self.add_var_opt('ephem-timecorr',et) self.__ephem_time = et def set_harmonics(self,h): # set model frequency harmonics self.add_var_opt('harmonics',h) self.__harmonics = h def set_Nlive(self,nl): # set number of live points self.add_var_opt('Nlive',nl) self.__Nlive = nl def set_Nmcmc(self,nm): # set number of MCMC iterations self.add_var_opt('Nmcmc',nm) self.__Nmcmc = nm def set_Nmcmcinitial(self,nm): # set number of MCMC iterations self.add_var_opt('Nmcmcinitial',nm) self.__Nmcmcinitial = nm def set_Nruns(self,nr): # set number of internal nested sample runs self.add_var_opt('Nruns',nr) self.__Nruns = nr def set_tolerance(self, tol): # set tolerance criterion for finishing nested sampling self.add_var_opt('tolerance', tol) self.__tolerance = tol def set_randomseed(self, rs): # set random number generator seed self.add_var_opt('randomseed', rs) self.__randomseed = rs def set_ensemble_stretch(self, f): # set the fraction of time to use ensemble stretch moves as proposal self.add_var_opt('ensembleStretch', f) self.__ensemble_stretch = f def set_ensemble_walk(self, f): # set the fraction of time to use ensemble walk moves as proposal self.add_var_opt('ensembleWalk', f) self.__ensemble_walk = f def set_temperature(self, temp): # set temperature scale for covariance proposal self.add_var_opt('temperature', temp) self.__temperature = temp def set_diffev(self,de): # set fraction of time to use differential evolution as proposal self.add_var_opt('diffev',de) self.__diffev = de def set_inject_file(self,ifil): # set a pulsar parameter file from which to make an injection self.add_var_opt('inject-file',ifil) self.__inject_file = ifil def set_inject_output(self,iout): # set filename in which to output the injected signal self.add_var_opt('inject-output',iout) self.__inject_output = iout def set_fake_data(self,fd): # set the detectors from which to generate fake data self.add_var_opt('fake-data',fd) self.__fake_data = fd def set_fake_psd(self,fp): # set the PSDs of the fake data self.add_var_opt('fake-psd',fp) self.__fake_psd = fp def set_fake_starts(self,fs): # set the start times of the fake data self.add_var_opt('fake-starts',fs) self.__fake_starts = fs def set_fake_lengths(self,fl): # set the lengths of the fake data self.add_var_opt('fake-lengths',fl) self.__fake_lengths = fl def set_fake_dt(self,fdt): # set the sample interval of the fake data self.add_var_opt('fake-dt',fdt) self.__fake_dt = fdt def set_scale_snr(self,ssnr): # set the SNR of the injected signal self.add_var_opt('scale-snr',ssnr) self.__scale_snr = ssnr def set_sample_files(self,ssf): # set the nested sample files to be used as a prior self.add_var_opt('sample-files',ssf) self.__sample_files = ssf def set_sample_nlives(self,snl): # set the number of live points for the nested sample files self.add_var_opt('sample-nlives',snl) self.__sample_nlives = snl def set_prior_cell(self,pc): # set the k-d tree cell size for the prior self.add_var_opt('prior-cell',pc) self.__prior_cell = pc def set_OldChunks(self): # use the old data segmentation routine i.e. 30 min segments self.add_var_opt('oldChunks', '') self.__oldChunks = True def set_source_model(self): # use the physical parameter model from Jones # add this into the generic 'macroargs' macro as it is not a value that is always required if 'macroargs' in self.get_opts(): curmacroval = self.get_opts()['macroargs'] else: curmacroval = '' curmacroval = curmacroval + ' --source-model' self.add_macro('macroargs', curmacroval) self.__sourceModel = True def set_biaxial(self): # the model is a biaxial star using the amplitude/phase waveform parameterisation # add this into the generic 'macroargs' macro as it is not a value that is always required if 'macroargs' in self.get_opts(): curmacroval = self.get_opts()['macroargs'] else: curmacroval = '' curmacroval = curmacroval + ' --biaxial' self.add_macro('macroargs', curmacroval) self.__biaxial = True def set_gaussian_like(self): # use a Gaussian likelihood rather than the Students-t likelihood self.add_var_opt('gaussian-like', '') self.__gaussian_like = True def set_randomise(self, f): # set flag to randomise the data times stamps for "background" analyses # add this into the generic 'macroargs' macro as it is not a value that is always required if 'macroargs' in self.get_opts(): curmacroval = self.get_opts()['macroargs'] else: curmacroval = '' curmacroval = curmacroval + ' --randomise ' + f self.add_macro('macroargs', curmacroval) self.__randomise = f def set_roq(self): # set to use Reduced Order Quadrature (ROQ) # add this into the generic 'macroargs' macro as it is not a value that is always required if 'macroargs' in self.get_opts(): curmacroval = self.get_opts()['macroargs'] else: curmacroval = '' curmacroval = curmacroval + ' --roq' self.add_macro('macroargs', curmacroval) self.__use_roq = True def set_roq_ntraining(self, f): # set the number of training waveforms to use in ROQ basis generation # add this into the generic 'macroargs' macro as it is not a value that is always required if 'macroargs' in self.get_opts(): curmacroval = self.get_opts()['macroargs'] else: curmacroval = '' curmacroval = curmacroval + ' --ntraining ' + f self.add_macro('macroargs', curmacroval) self.__roq_ntraining = f def set_roq_tolerance(self, f): # set the tolerance to use in ROQ basis generation # add this into the generic 'macroargs' macro as it is not a value that is always required if 'macroargs' in self.get_opts(): curmacroval = self.get_opts()['macroargs'] else: curmacroval = '' curmacroval = curmacroval + ' --roq-tolerance ' + f self.add_macro('macroargs', curmacroval) self.__roq_tolerance = f def set_roq_uniform(self): # set to use uniform distributions when sprinkling (phase) parameters for ROQ training basis generation # add this into the generic 'macroargs' macro as it is not a value that is always required if 'macroargs' in self.get_opts(): curmacroval = self.get_opts()['macroargs'] else: curmacroval = '' curmacroval = curmacroval + ' --roq-uniform' self.add_macro('macroargs', curmacroval) self.__roq_uniform = True def set_roq_inputweights(self,f): # set the location of the file containing pregenerated ROQ interpolants # add this into the generic 'macroargs' macro as it is not a value that is always required if 'macroargs' in self.get_opts(): curmacroval = self.get_opts()['macroargs'] else: curmacroval = '' curmacroval = curmacroval + ' --input-weights ' + f self.add_macro('macroargs', curmacroval) self.__roq_input_weights = f def set_roq_outputweights(self,f): # set the location of the file to output ROQ interpolants # add this into the generic 'macroargs' macro as it is not a value that is always required if 'macroargs' in self.get_opts(): curmacroval = self.get_opts()['macroargs'] else: curmacroval = '' curmacroval = curmacroval + ' --output-weights ' + f self.add_macro('macroargs', curmacroval) self.__roq_output_weights = f def set_roq_chunkmax(self,f): # set the maximum chunk length for if using ROQ (just adding using the chunk-max argument) # add this into the generic 'macroargs' macro as it is not a value that is always required if 'macroargs' in self.get_opts(): curmacroval = self.get_opts()['macroargs'] else: curmacroval = '' curmacroval = curmacroval + ' --chunk-max ' + f self.add_macro('macroargs', curmacroval) self.__chunk_max = int(f) """ Job for creating the result page for a particular source """ class resultpageJob(pipeline.CondorDAGJob, pipeline.AnalysisJob): def __init__(self, execu, univ='local', accgroup=None, accuser=None, logdir=None, rundir=None): self.__executable = execu self.__universe = univ pipeline.CondorDAGJob.__init__(self, self.__universe, self.__executable) pipeline.AnalysisJob.__init__(self, None) if accgroup != None: self.add_condor_cmd('accounting_group', accgroup) if accuser != None: self.add_condor_cmd('accounting_group_user', accuser) self.add_condor_cmd('getenv','True') # set log files for job if logdir != None: self.set_stdout_file(os.path.join(logdir, 'resultpage-$(cluster).out')) self.set_stderr_file(os.path.join(logdir, 'resultpage-$(cluster).err')) else: self.set_stdout_file('resultpage-$(cluster).out') self.set_stderr_file('resultpage-$(cluster).err') if rundir != None: self.set_sub_file(os.path.join(rundir, 'resultpage.sub')) else: self.set_sub_file('resultpage.sub') self.add_arg('$(macroconfigfile)') # macro for input configuration file class resultpageNode(pipeline.CondorDAGNode, pipeline.AnalysisNode): """ A resultpageNode runs an instance of the result page script in a condor DAG. """ def __init__(self,job): """ job = A CondorDAGJob that can run an instance of lalapps_knope_result_page """ pipeline.CondorDAGNode.__init__(self,job) pipeline.AnalysisNode.__init__(self) self.__configfile = None def set_config(self, configfile): self.add_macro('macroconfigfile', configfile) self.__configfile = configfile """ Job for creating the collated results page for all sources """ class collateJob(pipeline.CondorDAGJob, pipeline.AnalysisJob): def __init__(self, execu, univ='local', accgroup=None, accuser=None, logdir=None, rundir=None): self.__executable = execu self.__universe = univ pipeline.CondorDAGJob.__init__(self, self.__universe, self.__executable) pipeline.AnalysisJob.__init__(self, None) if accgroup != None: self.add_condor_cmd('accounting_group', accgroup) if accuser != None: self.add_condor_cmd('accounting_group_user', accuser) self.add_condor_cmd('getenv','True') # set log files for job if logdir != None: self.set_stdout_file(os.path.join(logdir, 'collate-$(cluster).out')) self.set_stderr_file(os.path.join(logdir, 'collate-$(cluster).err')) else: self.set_stdout_file('collate-$(cluster).out') self.set_stderr_file('collate-$(cluster).err') if rundir != None: self.set_sub_file(os.path.join(rundir, 'collate.sub')) else: self.set_sub_file('collate.sub') self.add_arg('$(macroconfigfile)') # macro for input configuration file class collateNode(pipeline.CondorDAGNode, pipeline.AnalysisNode): """ A collateNode runs an instance of the result page collation script in a condor DAG. """ def __init__(self,job): """ job = A CondorDAGJob that can run an instance of lalapps_knope_collate_results """ pipeline.CondorDAGNode.__init__(self,job) pipeline.AnalysisNode.__init__(self) self.__configfile = None def set_config(self, configfile): self.add_macro('macroconfigfile', configfile) self.__configfile = configfile class nest2posJob(pipeline.CondorDAGJob, pipeline.AnalysisJob): """ A merge nested sampling files job to use lalinference_nest2pos """ def __init__(self, execu, univ='local', accgroup=None, accuser=None, logdir=None, rundir=None): self.__executable = execu self.__universe = univ pipeline.CondorDAGJob.__init__(self, self.__universe, self.__executable) pipeline.AnalysisJob.__init__(self, None) if accgroup != None: self.add_condor_cmd('accounting_group', accgroup) if accuser != None: self.add_condor_cmd('accounting_group_user', accuser) self.add_condor_cmd('getenv','True') # set log files for job if logdir != None: self.set_stdout_file(os.path.join(logdir, 'n2p-$(cluster).out')) self.set_stderr_file(os.path.join(logdir, 'n2p-$(cluster).err')) else: self.set_stdout_file('n2p-$(cluster).out') self.set_stderr_file('n2p-$(cluster).err') self.add_arg('--non-strict-versions') # force use of --non-strict-versions flag self.add_arg('$(macroinputfiles)') # macro for input nested sample files if rundir != None: self.set_sub_file(os.path.join(rundir, 'n2p.sub')) else: self.set_sub_file('n2p.sub') class nest2posNode(pipeline.CondorDAGNode, pipeline.AnalysisNode): """ A nest2posNode runs a instance of the lalinference_nest2pos to combine individual nested sample files in a condor DAG. """ def __init__(self,job): """ job = A CondorDAGJob that can run an of the nested sample combination script """ pipeline.CondorDAGNode.__init__(self,job) pipeline.AnalysisNode.__init__(self) self.set_retry(1) # retry the node once # initilise job variables self.__nest_files = None self.__nest_live = None self.__outfile = None self.__header = None self.__npos = None self.__gzip = False def set_nest_files(self,nestfiles): # set all the nested sample files self.__nest_files = nestfiles fe = os.path.splitext(nestfiles[0])[-1].lower() # set header file (only if not using hdf5 output) if fe != '.hdf' and fe != '.h5': header = nestfiles[0].rstrip('.gz')+'_params.txt' self.__header = header self.add_var_opt('headers', header) self.add_macro('macroinputfiles', ' '.join(nestfiles)) def set_nest_live(self,nestlive): # set the number of live points from each file self.add_var_opt('Nlive',nestlive) self.__nest_live = nestlive def set_outfile(self,outfile): # set the output file for the posterior file self.add_var_opt('pos', outfile) self.__outfile = outfile def set_numpos(self, npos): # set the number of posterior samples for the posterior file self.add_var_opt('npos', npos) self.__npos = npos def set_gzip(self): self.add_var_opt('gzip', '') self.__gzip = True # DEPRECATED class createresultspageJob(pipeline.CondorDAGJob, pipeline.AnalysisJob): """ A job to create an individual pulsar results page """ def __init__(self,execu,logpath,accgroup,accuser): self.__executable = execu self.__universe = 'vanilla' pipeline.CondorDAGJob.__init__(self, self.__universe, self.__executable) pipeline.AnalysisJob.__init__(self, None) self.add_condor_cmd('getenv','True') self.add_condor_cmd('accounting_group', accgroup) self.add_condor_cmd('accounting_group_user', accuser) self.set_stdout_file(logpath+'/create_results_page-$(cluster).out') self.set_stderr_file(logpath+'/create_results_page-$(cluster).err') self.set_sub_file('create_results_page.sub') # additional required args self.add_arg('$(macrom)') # macro for MCMC directories self.add_arg('$(macrobk)') # macro for Bk (fine heterodyne file) directories self.add_arg('$(macroi)') # macro for IFOs self.add_arg('$(macrof)') # macro for nested sampling files self.add_arg('$(macrow)') # macro to say if a hardware injection self.add_arg('$(macros)') # macro to say if a software injection class createresultspageNode(pipeline.CondorDAGNode, pipeline.AnalysisNode): """ A createresultspage node to run as part of a condor DAG. """ def __init__(self,job): """ job = A CondorDAGJob that can run the segment list finding script """ pipeline.CondorDAGNode.__init__(self,job) pipeline.AnalysisNode.__init__(self) # initilise job variables self.__outpath = None self.__domcmc = False self.__mcmcdirs = [] self.__donested = False self.__nestedfiles = [] self.__parfile = None self.__Bkfiles = [] self.__priorfile = None self.__ifos = [] self.__histbins = None self.__epsout = False def set_outpath(self,val): # set the detector self.add_var_opt('o', val, short=True) self.__outpath = val def set_domcmc(self): # set to say using MCMC chains as input self.add_var_opt('M', '', short=True) self.__domcmc = True def set_mcmcdir(self,val): # set the MCMC file directories macroval = '' for f in val: macroval = '%s-m %s ' % (macroval, f) self.add_macro('macrom', macroval) self.add_macro('macrof', '') # empty macro for nested files self.__mcmcdirs = val def set_donested(self): # set to say using nested sampling results as input self.add_var_opt('nested', '') self.__donested = True def set_nestedfiles(self,val): # set the nested sampling files macroval = '' for f in val: macroval = '%s-f %s ' % (macroval, f) self.add_macro('macrof', macroval) self.add_macro('macrom', '') # empty macro for mcmc directories self.__nestedfiles = val def set_parfile(self,val): # set the pulsar parameter file self.add_var_opt('p', val, short=True) self.__parfile = val def set_bkfiles(self,val): # set the fine heterodyned data files macroval = '' for f in val: macroval = '%s-b %s ' % (macroval, f) self.add_macro('macrobk', macroval) self.__Bkfiles = val def set_priordir(self,val): # set the prior file directory self.add_var_opt('r', val, short=True) self.__priordir = None def set_ifos(self,val): # set the IFOs to analyse macroval = '' for f in val: macroval = '%s-i %s ' % (macroval, f) self.add_macro('macroi', macroval) self.__ifos = val def set_histbins(self,val): # set the number of histogram bins self.add_var_opt('n', val, short=True) self.__histbins = val def set_epsout(self): # set to output eps figs self.add_var_opt('e', '', short=True) self.__epsout = True def set_swinj(self, isswinj): # set to say that analysing software injection if isswinj: self.add_macro('macros', '--sw-inj') else: self.add_macro('macros', '') def set_hwinj(self, ishwinj): # set to say that analysing hardware injection if ishwinj: self.add_macro('macrow', '--hw-inj') else: self.add_macro('macrow', '') class collateresultsJob(pipeline.CondorDAGJob, pipeline.AnalysisJob): """ A job to collate all the individual pulsar results pages """ def __init__(self,execu,logpath,accgroup,accuser): self.__executable = execu self.__universe = 'vanilla' pipeline.CondorDAGJob.__init__(self, self.__universe, self.__executable) pipeline.AnalysisJob.__init__(self, None) self.add_condor_cmd('getenv','True') self.add_condor_cmd('accounting_group', accgroup) self.add_condor_cmd('accounting_group_user', accuser) self.set_stdout_file(logpath+'/collate_results-$(cluster).out') self.set_stderr_file(logpath+'/collate_results-$(cluster).err') self.set_sub_file('collate_results.sub') # some required argument macros self.add_arg('$(macroifo)') # for IFOs #self.add_arg('$(macrou)') # for output upper limits #self.add_arg('$(macron)') # for output pulsar values class collateresultsNode(pipeline.CondorDAGNode, pipeline.AnalysisNode): """ A collateresults node to run as part of a condor DAG. """ def __init__(self,job): """ job = A CondorDAGJob that can run the segment list finding script """ pipeline.CondorDAGNode.__init__(self,job) pipeline.AnalysisNode.__init__(self) # initilise job variables self.__outpath = None self.__inpath = None self.__parfile = None self.__compilelatex = False self.__sorttype = None self.__ifos = [] self.__outputlims = [] self.__outputvals = [] self.__outputhist = False self.__outputulplot = False self.__withprior = False self.__epsout = False def set_outpath(self,val): # set the detector self.add_var_opt('o', val, short=True) self.__outpath = val def set_inpath(self,val): # set the input path self.add_var_opt('z', val, short=True) self.__inpath = val def set_parfile(self,val): # set the pulsar parameter file directory self.add_var_opt('p', val, short=True) self.__parfile = val def set_compilelatex(self): # set to compile LaTeX results table self.add_var_opt('l', '', short=True) self.__compilelatex = True def set_sorttype(self,val): # set the sorting order of the output results self.add_var_opt('s', val, short=True) self.__sorttype = val def set_ifos(self,val): # set the list if IFOs to output results for macroval = '' for f in val: macroval = '%s-i %s ' % (macroval, f) self.add_macro('macroifo', macroval) self.__ifos = val def set_outputlims(self,val): # set the upper limit results to output macroval = '' for f in val: macroval = '%s-u %s ' % (macroval, f) self.add_macro('macrou', macroval) self.__outputlims = val def set_outputvals(self,val): # set the pulsar parameter values to output macroval = '' for f in val: macroval = '%s-n %s ' % (macroval, f) self.add_macro('macron', macroval) self.__outputvals = val def set_outputhist(self): # set to output histograms of the results self.add_var_opt('k', '', short=True) self.__outputhist = True def set_outputulplot(self): # set to output a plot of the ULs self.add_var_opt('t', '', short=True) self.__outputulplot = True def set_withprior(self): # set to output prior values with the hist and UL plots self.add_var_opt('w', '', short=True) self.__withprior = True def set_epsout(self): # set to output plots in eps format as well as png self.add_var_opt('e', '', short=True) self.__epsout = True