Loading minke/mdctools.py +84 −30 Original line number Diff line number Diff line Loading @@ -7,6 +7,7 @@ from pylal.antenna import response from pylal.date import XLALTimeDelayFromEarthCenter from pylal.xlal.datatypes.ligotimegps import LIGOTimeGPS from pylal import inject from glue.ligolw.utils import process import lal, lalframe from pylal import Fr Loading Loading @@ -53,41 +54,25 @@ class MDCSet(): waveforms = [] def __init__(self, detectors, simtable, name=None, full=True): def __init__(self, detectors, name='MDC Set'): """ Represents an MDC set, stored in an XML SimBurstTable file. Parameters ---------- detectors : list A list of detector names where the injections should be made simtable : str The filepath to a simbursttable xml file. A list of detector names where the injections should be made. name : str A name for the MDC Set. Defaults to 'MDC Set'. """ self.detectors = detectors sim_burst_table = lalburst.SimBurstTableFromLIGOLw(simtable, None, None) self.waveforms = [] self.strains = [] self.egw = [] self.times = [] while True: # This is an ugly kludge to get around the poor choice of wavform name in the xmls, and if sim_burst_table.waveform[:3]=="s15": self.numrel_file = str(sim_burst_table.waveform) sim_burst_table.waveform = "Dimmelmeier+08" self.waveforms.append(sim_burst_table) if full: self._generate_burst(sim_burst_table) self._measure_hrss() self._measure_egw_rsq() self.times.append(sim_burst_table.time_geocent_gps) if sim_burst_table.next is None: break sim_burst_table = sim_burst_table.next if name: self.name = name else: self.name = self._simID(0) self.times = np.array(self.times) Loading @@ -101,11 +86,77 @@ class MDCSet(): The waveform which should be added to the MDC set. """ self.waveforms += waveform self.times += waveform.time_geocent_gps() self.waveforms.append(waveform) self.times = np.append(self.times, waveform.time) def save_xml(self): pass def save_xml(self, filename): """ Save the MDC set as an XML SimBurstTable. Parameters ---------- filename : str The location to save the xml file. The output is gzipped, so ending it with a ".gz" would stick with convention. """ xmldoc = ligolw.Document() lw = xmldoc.appendChild(ligolw.LIGO_LW()) sim = lsctables.New(lsctables.SimBurstTable) lw.appendChild(sim) # This needs to be given the proper metadata once the package has the maturity to # write something sensible. procrow = process.register_to_xmldoc(xmldoc, "pyburst_binj", {}) for waveform in self.waveforms: sim.append(waveform._row(sim)) # Write out the xml and gzip it. utils.write_filename(xmldoc, filename, gz=True) def load_xml(self, filename, full=True, start=None, stop=None): """Load the MDC Set from an XML file containing the SimBurstTable. Parameters ---------- filename : str The filename of the XML file. full : bool If this is true (which is the default) then all of the calculated parameters are computed from the waveform definintion. start : float The time at which the xml read-in should start. The default is "None", in which case the xml file will be read-in from the start. end : float The last time to be read from the xml file. The default is None, which causes the xml to be read right-up to the last time in the file. To Do ----- A the moment this loads the information in to the object, but it doesn't produce waveform objects for each of the injections in the file. This should be fixed so that the object works symmetrically. """ sim_burst_table = lalburst.SimBurstTableFromLIGOLw(filename, start, stop) while True: # This is an ugly kludge to get around the poor choice of wavform name in the xmls, and if sim_burst_table.waveform[:3]=="s15": self.numrel_file = str(sim_burst_table.waveform) sim_burst_table.waveform = "Dimmelmeier+08" self.waveforms.append(sim_burst_table) if full: self._generate_burst(sim_burst_table) self._measure_hrss() self._measure_egw_rsq() self.times.append(sim_burst_table.time_geocent_gps) if sim_burst_table.next is None: break sim_burst_table = sim_burst_table.next del(sim_burst_table) def _generate_burst(self, row,rate=16384.0): """ Loading @@ -114,8 +165,8 @@ class MDCSet(): Parameters ---------- row : int The row number of the waveforms to be measured row : SimBurst Row The row of the waveform to be measured rate : float The sampling rate of the signal, in Hz. Defaults to 16384.0Hz Loading @@ -140,7 +191,10 @@ class MDCSet(): print a, getattr(row,a) continue # the structure is different than the TableRow #print row.numrel_data try: self.swig_row.numrel_data = row.numrel_data except: pass hp, hx = lalburst.GenerateSimBurst(self.swig_row, 1.0/rate) # FIXME: Totally inefficent --- but can we deep copy a SWIG SimBurst? # DW: I tried that, and it doesn't seem to work :/ Loading minke/sources.py +74 −7 Original line number Diff line number Diff line Loading @@ -20,6 +20,7 @@ import lalsimulation from minke.distribution import * import matplotlib.pyplot as plt class Waveform(object): """ Loading Loading @@ -108,6 +109,15 @@ class Waveform(object): return pol_ellipse_e, pol_ellipse_angle def plot(self): """ Produce a plot of the injection. """ self._generate() plt.plot(self.hp.data.data, label="+ polarisation") plt.plot(self.hx.data.data, label="x polarisation") def uniform_time(self): """ Get a set of randomized (integer) event times. Loading @@ -121,18 +131,22 @@ class Waveform(object): Todo ---- This can currently only make injections on a uniform sky. This should be fixed to take a generic distribution function. We also need to add the process_id back in. """ self.row = sim.RowType() # Required columns not defined makes ligolw unhappy for a in lsctables.SimBurstTable.validcolumns.keys(): setattr(self.row, a, None) self.row.waveform = self.waveform self.row.set_time_geocent(self.time) self.row.set_time_geocent(GPS(float(self.time))) # Right now this only does uniform sky distributions, but we should provide a way to do /any/ distribution. self.row.ra, self.row.dec, self.row.psi = self.uniform_sky() self.row.simulation_id = sim.get_next_id() self.row.waveform_number = random.randint(0,int(2**32)-1) self.row.process_id = procrow.process_id ### !! This needs to be updated. self.row.process_id = "process:process_id:0" #procrow.process_id self.row.time_slide_id = ilwd.ilwdchar("time_slide:time_slide_id:%d" % options.time_slide_id) Loading Loading @@ -199,10 +213,14 @@ class SineGaussian(Waveform): self.tstart, self.tstop = time[0], time[1] self.seed = seed random.seed(seed) self.time = random.randint(self.tstart, self.tstop, 1) + random.rand(1) self.time = (random.randint(self.tstart, self.tstop, 1) + random.rand(1))[0] else: self.time = time self.polarisation = polarisation self.q, self.frequency, self.hrss = self.draw_params() def draw_params(self): """ Draw a set of parameters (hrss, q, frequency) appropriate for a set of sinegaussian type injections. Loading Loading @@ -237,26 +255,75 @@ class SineGaussian(Waveform): return q, f, h def _row(self, sim): def _generate(self, rate=16384.0): """ Generate the burst described in a given row, so that it can be measured. Parameters ---------- rate : float The sampling rate of the signal, in Hz. Defaults to 16384.0Hz Returns ------- hp : The strain in the + polarisation hx : The strain in the x polarisation hp0 : A copy of the strain in the + polarisation hx0 : A copy of the strain in the x polarisation """ row = self._row() self.swig_row = lalburst.CreateSimBurst() for a in lsctables.SimBurstTable.validcolumns.keys(): try: setattr(self.swig_row, a, getattr( row, a )) except AttributeError: continue # we didn't define it except TypeError: #print a, getattr(row,a) continue # the structure is different than the TableRow #print row.numrel_data try: self.swig_row.numrel_data = row.numrel_data except: pass hp, hx = lalburst.GenerateSimBurst(self.swig_row, 1.0/rate) # FIXME: Totally inefficent --- but can we deep copy a SWIG SimBurst? # DW: I tried that, and it doesn't seem to work :/ hp0, hx0 = lalburst.GenerateSimBurst(self.swig_row, 1.0/rate) self.hp, self.hx, self.hp0, self.hx0 = hp, hx, hp0, hx0 def _row(self, sim=None): """ Produce a simburst table row for this waveform. Parameters ---------- sim : lsctable The sim table for this row. Todo ---- This can currently only make injections on a uniform sky. This should be fixed to take a generic distribution function. Need to make sure that the correct process ID is used. """ if not sim: sim = lsctables.New(lsctables.SimBurstTable) self.row = sim.RowType() # Required columns not defined makes ligolw unhappy for a in lsctables.SimBurstTable.validcolumns.keys(): setattr(self.row, a, None) self.row.waveform = self.waveform self.row.set_time_geocent(self.time) self.row.set_time_geocent(GPS(float(self.time))) # Right now this only does uniform sky distributions, but we should provide a way to do /any/ distribution. self.row.ra, self.row.dec, self.row.psi = self.uniform_sky() self.row.simulation_id = sim.get_next_id() self.row.waveform_number = random.randint(0,int(2**32)-1) self.row.process_id = procrow.process_id self.row.time_slide_id = ilwd.ilwdchar("time_slide:time_slide_id:%d" % options.time_slide_id) self.row.process_id = "process:process_id:0" #procrow.process_id self.row.time_slide_id = ilwd.ilwdchar("time_slide:time_slide_id:%d" % 1)#options.time_slide_id) self.row.q, self.row.frequency, self.row.hrss = self.draw_params() self.row.pol_ellipse_e, self.row.pol_ellipse_angle = self.parse_polarisation(self.polarisation) Loading Loading
minke/mdctools.py +84 −30 Original line number Diff line number Diff line Loading @@ -7,6 +7,7 @@ from pylal.antenna import response from pylal.date import XLALTimeDelayFromEarthCenter from pylal.xlal.datatypes.ligotimegps import LIGOTimeGPS from pylal import inject from glue.ligolw.utils import process import lal, lalframe from pylal import Fr Loading Loading @@ -53,41 +54,25 @@ class MDCSet(): waveforms = [] def __init__(self, detectors, simtable, name=None, full=True): def __init__(self, detectors, name='MDC Set'): """ Represents an MDC set, stored in an XML SimBurstTable file. Parameters ---------- detectors : list A list of detector names where the injections should be made simtable : str The filepath to a simbursttable xml file. A list of detector names where the injections should be made. name : str A name for the MDC Set. Defaults to 'MDC Set'. """ self.detectors = detectors sim_burst_table = lalburst.SimBurstTableFromLIGOLw(simtable, None, None) self.waveforms = [] self.strains = [] self.egw = [] self.times = [] while True: # This is an ugly kludge to get around the poor choice of wavform name in the xmls, and if sim_burst_table.waveform[:3]=="s15": self.numrel_file = str(sim_burst_table.waveform) sim_burst_table.waveform = "Dimmelmeier+08" self.waveforms.append(sim_burst_table) if full: self._generate_burst(sim_burst_table) self._measure_hrss() self._measure_egw_rsq() self.times.append(sim_burst_table.time_geocent_gps) if sim_burst_table.next is None: break sim_burst_table = sim_burst_table.next if name: self.name = name else: self.name = self._simID(0) self.times = np.array(self.times) Loading @@ -101,11 +86,77 @@ class MDCSet(): The waveform which should be added to the MDC set. """ self.waveforms += waveform self.times += waveform.time_geocent_gps() self.waveforms.append(waveform) self.times = np.append(self.times, waveform.time) def save_xml(self): pass def save_xml(self, filename): """ Save the MDC set as an XML SimBurstTable. Parameters ---------- filename : str The location to save the xml file. The output is gzipped, so ending it with a ".gz" would stick with convention. """ xmldoc = ligolw.Document() lw = xmldoc.appendChild(ligolw.LIGO_LW()) sim = lsctables.New(lsctables.SimBurstTable) lw.appendChild(sim) # This needs to be given the proper metadata once the package has the maturity to # write something sensible. procrow = process.register_to_xmldoc(xmldoc, "pyburst_binj", {}) for waveform in self.waveforms: sim.append(waveform._row(sim)) # Write out the xml and gzip it. utils.write_filename(xmldoc, filename, gz=True) def load_xml(self, filename, full=True, start=None, stop=None): """Load the MDC Set from an XML file containing the SimBurstTable. Parameters ---------- filename : str The filename of the XML file. full : bool If this is true (which is the default) then all of the calculated parameters are computed from the waveform definintion. start : float The time at which the xml read-in should start. The default is "None", in which case the xml file will be read-in from the start. end : float The last time to be read from the xml file. The default is None, which causes the xml to be read right-up to the last time in the file. To Do ----- A the moment this loads the information in to the object, but it doesn't produce waveform objects for each of the injections in the file. This should be fixed so that the object works symmetrically. """ sim_burst_table = lalburst.SimBurstTableFromLIGOLw(filename, start, stop) while True: # This is an ugly kludge to get around the poor choice of wavform name in the xmls, and if sim_burst_table.waveform[:3]=="s15": self.numrel_file = str(sim_burst_table.waveform) sim_burst_table.waveform = "Dimmelmeier+08" self.waveforms.append(sim_burst_table) if full: self._generate_burst(sim_burst_table) self._measure_hrss() self._measure_egw_rsq() self.times.append(sim_burst_table.time_geocent_gps) if sim_burst_table.next is None: break sim_burst_table = sim_burst_table.next del(sim_burst_table) def _generate_burst(self, row,rate=16384.0): """ Loading @@ -114,8 +165,8 @@ class MDCSet(): Parameters ---------- row : int The row number of the waveforms to be measured row : SimBurst Row The row of the waveform to be measured rate : float The sampling rate of the signal, in Hz. Defaults to 16384.0Hz Loading @@ -140,7 +191,10 @@ class MDCSet(): print a, getattr(row,a) continue # the structure is different than the TableRow #print row.numrel_data try: self.swig_row.numrel_data = row.numrel_data except: pass hp, hx = lalburst.GenerateSimBurst(self.swig_row, 1.0/rate) # FIXME: Totally inefficent --- but can we deep copy a SWIG SimBurst? # DW: I tried that, and it doesn't seem to work :/ Loading
minke/sources.py +74 −7 Original line number Diff line number Diff line Loading @@ -20,6 +20,7 @@ import lalsimulation from minke.distribution import * import matplotlib.pyplot as plt class Waveform(object): """ Loading Loading @@ -108,6 +109,15 @@ class Waveform(object): return pol_ellipse_e, pol_ellipse_angle def plot(self): """ Produce a plot of the injection. """ self._generate() plt.plot(self.hp.data.data, label="+ polarisation") plt.plot(self.hx.data.data, label="x polarisation") def uniform_time(self): """ Get a set of randomized (integer) event times. Loading @@ -121,18 +131,22 @@ class Waveform(object): Todo ---- This can currently only make injections on a uniform sky. This should be fixed to take a generic distribution function. We also need to add the process_id back in. """ self.row = sim.RowType() # Required columns not defined makes ligolw unhappy for a in lsctables.SimBurstTable.validcolumns.keys(): setattr(self.row, a, None) self.row.waveform = self.waveform self.row.set_time_geocent(self.time) self.row.set_time_geocent(GPS(float(self.time))) # Right now this only does uniform sky distributions, but we should provide a way to do /any/ distribution. self.row.ra, self.row.dec, self.row.psi = self.uniform_sky() self.row.simulation_id = sim.get_next_id() self.row.waveform_number = random.randint(0,int(2**32)-1) self.row.process_id = procrow.process_id ### !! This needs to be updated. self.row.process_id = "process:process_id:0" #procrow.process_id self.row.time_slide_id = ilwd.ilwdchar("time_slide:time_slide_id:%d" % options.time_slide_id) Loading Loading @@ -199,10 +213,14 @@ class SineGaussian(Waveform): self.tstart, self.tstop = time[0], time[1] self.seed = seed random.seed(seed) self.time = random.randint(self.tstart, self.tstop, 1) + random.rand(1) self.time = (random.randint(self.tstart, self.tstop, 1) + random.rand(1))[0] else: self.time = time self.polarisation = polarisation self.q, self.frequency, self.hrss = self.draw_params() def draw_params(self): """ Draw a set of parameters (hrss, q, frequency) appropriate for a set of sinegaussian type injections. Loading Loading @@ -237,26 +255,75 @@ class SineGaussian(Waveform): return q, f, h def _row(self, sim): def _generate(self, rate=16384.0): """ Generate the burst described in a given row, so that it can be measured. Parameters ---------- rate : float The sampling rate of the signal, in Hz. Defaults to 16384.0Hz Returns ------- hp : The strain in the + polarisation hx : The strain in the x polarisation hp0 : A copy of the strain in the + polarisation hx0 : A copy of the strain in the x polarisation """ row = self._row() self.swig_row = lalburst.CreateSimBurst() for a in lsctables.SimBurstTable.validcolumns.keys(): try: setattr(self.swig_row, a, getattr( row, a )) except AttributeError: continue # we didn't define it except TypeError: #print a, getattr(row,a) continue # the structure is different than the TableRow #print row.numrel_data try: self.swig_row.numrel_data = row.numrel_data except: pass hp, hx = lalburst.GenerateSimBurst(self.swig_row, 1.0/rate) # FIXME: Totally inefficent --- but can we deep copy a SWIG SimBurst? # DW: I tried that, and it doesn't seem to work :/ hp0, hx0 = lalburst.GenerateSimBurst(self.swig_row, 1.0/rate) self.hp, self.hx, self.hp0, self.hx0 = hp, hx, hp0, hx0 def _row(self, sim=None): """ Produce a simburst table row for this waveform. Parameters ---------- sim : lsctable The sim table for this row. Todo ---- This can currently only make injections on a uniform sky. This should be fixed to take a generic distribution function. Need to make sure that the correct process ID is used. """ if not sim: sim = lsctables.New(lsctables.SimBurstTable) self.row = sim.RowType() # Required columns not defined makes ligolw unhappy for a in lsctables.SimBurstTable.validcolumns.keys(): setattr(self.row, a, None) self.row.waveform = self.waveform self.row.set_time_geocent(self.time) self.row.set_time_geocent(GPS(float(self.time))) # Right now this only does uniform sky distributions, but we should provide a way to do /any/ distribution. self.row.ra, self.row.dec, self.row.psi = self.uniform_sky() self.row.simulation_id = sim.get_next_id() self.row.waveform_number = random.randint(0,int(2**32)-1) self.row.process_id = procrow.process_id self.row.time_slide_id = ilwd.ilwdchar("time_slide:time_slide_id:%d" % options.time_slide_id) self.row.process_id = "process:process_id:0" #procrow.process_id self.row.time_slide_id = ilwd.ilwdchar("time_slide:time_slide_id:%d" % 1)#options.time_slide_id) self.row.q, self.row.frequency, self.row.hrss = self.draw_params() self.row.pol_ellipse_e, self.row.pol_ellipse_angle = self.parse_polarisation(self.polarisation) Loading
