Added support for Gaussians, and converging very fast to the complete... (4fcf46e9) · Commits · Daniel Williams / Minke

.travis.yml

+1 −0

Original line number	Diff line number	Diff line
		@@ -15,6 +15,7 @@ addons:
		- libgsl0-dev
		- swig
		- bc
		- fftw-dev
		# nds2 dependencies
		- libsasl2-2
		# misc python dependencies

minke/distribution.py

+102 −1

Original line number	Diff line number	Diff line
		@@ -31,10 +31,111 @@ def uniform_phi(num):
		"""
		return random.random_sample(num) * 2 * numpy.pi

		def uniform_sky(number=1):
		"""
		Get a set of (RA, declination, polarization) randomized appopriately to astrophysical sources isotropically distributed in the sky.
		"""
		expnum = number
		ra = uniform_phi(expnum)
		dec = uniform_dec(expnum)
		pol = uniform_phi(expnum)
		return ra, dec, pol

		def favorable_sky(net, time):
		"""
		Wander through the skies, searching for a most favorable location --- draw extrinsic parameters as if the network antenna pattern magnitude were the PDF.
		"""
		ndraws = len(time)
		ra_out, dec_out, psi_out = numpy.empty((3, len(time)))
		while ndraws > 0:
		ra = numpy.random.uniform(0, 2 * numpy.pi, 1000)
		dec = numpy.random.uniform(-numpy.pi / 2, numpy.pi / 2, 1000)
		psi = numpy.random.uniform(0, 2 * numpy.pi, 1000)
		rnd = numpy.random.uniform(0, 1, 1000)
		for r, d, p, n in zip(ra, dec, psi, rnd):
		if n < sky_params(net, time, r, d, p):
		ndraws -= 1
		ra_out[ndraws] = r
		dec_out[ndraws] = d
		psi_out[ndraws] = p
		break
		return ra_out, dec_out, psi_out

		def uniform_interval(interval, num):
		"""
		Uniform in an interval, with the interval being a 2-tuple. num controls the number of draws. See also:
		Return a number, or a list of numbers which are sampled
		from a uniform distribution.

		Parameters
		----------
		interval : tuple (lower, upper)
		The interval which the uniform distribution covers.
		num : int
		The number of samples which should be drawn from the distribution.

		Returns
		-------
		sample : float or list
		A sample, or a list of samples.

		Notes
		-----

		http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.random_sample.html#numpy.random.random_sample
		"""
		return (interval[1] - interval[0]) * random.random_sample(num) + interval[0]

		def uniform_time(tstart, tstop, number):
		"""
		Get a set of randomized (integer) event times.
		"""
		return random.randint(tstart, tstop, number) + random.rand(number)

		def log_uniform(lower, upper, number):
		"""
		Draw uniformly in the log of a predefined range.

		Parameters
		----------
		lower : float
		The lower hrss value (n.b. not the lower log(hrss) value)
		upper : float
		The upper hrss value (n.b. not the upper log(hrss) value)
		number : int
		The number of samples to be drawn from the distribution.

		Returns
		-------
		sample : float
		An sample value drawn from the log uniform distribution.
		"""
		log10h = (numpy.log10(lower), numpy.log10(upper))
		return 10**uniform_interval(log10h, number)

		def even_time(start, stop, rate, jitter=0):
		"""
		Produce an evenly-distributed set of times which has some jitter.

		Parameters
		----------
		start : int
		The gps start time of the set.
		stop : int
		The gps end time of the set.
		rate : int
		The rate of events per year.
		jitter : float
		The number of seconds around which the time can "jitter".
		This has the effect of adding a small random time on to each
		time returned in the distribution, up to a maximum of half
		the jitter value, or taking away up to half the jitter value.

		Returns
		-------
		times : ndarray
		A numpy array of gps times where injections should be made.
		"""
		expnum_exact = (stop-start) * rate / 365.0 / 24 / 3600
		interval = (stop - start) / expnum_exact
		time = numpy.array(map(lambda i: start + iinterval + jitter (random.rand()-0.5), range(1,int(expnum_exact)+1)))
		return time

minke/mdctools.py

+3 −1

Original line number	Diff line number	Diff line
		@@ -152,7 +152,9 @@ class MDCSet():
		self._generate_burst(sim_burst_table)
		self._measure_hrss()
		self._measure_egw_rsq()
		self.times.append(sim_burst_table.time_geocent_gps)

		self.times = np.append(self.times, float(sim_burst_table.time_geocent_gps))
		#np.insert(self.times, len(self.times), sim_burst_table.time_geocent_gps)
		if sim_burst_table.next is None: break
		sim_burst_table = sim_burst_table.next

minke/sources.py

+110 −200

Original line number	Diff line number	Diff line
		@@ -29,53 +29,17 @@ class Waveform(object):
		In the future, different sources should subclass this and override the generation routines.
		"""

		sim = lsctables.New(lsctables.SimBurstTable)

		numrel_data = []
		waveform = "Generic"
		expnum = 1

		def log_hrss(self):
		"""
		Draw uniformly in the log of a predefined hrss range
		"""
		log10h = segment(numpy.log10(self.h_values[0]), numpy.log10(self.h_values[1]))
		return 10**uniform_interval(log10h, 1)

		def log_egw(self):
		"""
		Draw uniformly in the log of a predefined E_gw range
		"""
		log10e = segment(numpy.log10(self.egw_range[0]), numpy.log10(self.egw_range[1]))
		return 10**uniform_interval(log10e, self.expnum)

		def uniform_sky(self):
		"""
		Get a set of (RA, declination, polarization) randomized appopriately to astrophysical sources isotropically distributed in the sky.
		"""
		expnum = self.expnum
		ra = uniform_phi(expnum)
		dec = uniform_dec(expnum)
		pol = uniform_phi(expnum)
		return ra, dec, pol
		def _clear_params(self):
		self.params = {}
		for a in lsctables.SimBurstTable.validcolumns.keys():
		self.params[a] = None

		def favorable_sky(net, time):
		"""
		Wander through the skies, searching for a most favorable location --- draw extrinsic parameters as if the network antenna pattern magnitude were the PDF.
		"""
		ndraws = len(time)
		ra_out, dec_out, psi_out = numpy.empty((3, len(time)))
		while ndraws > 0:
		ra = numpy.random.uniform(0, 2 * numpy.pi, 1000)
		dec = numpy.random.uniform(-numpy.pi / 2, numpy.pi / 2, 1000)
		psi = numpy.random.uniform(0, 2 * numpy.pi, 1000)
		rnd = numpy.random.uniform(0, 1, 1000)
		for r, d, p, n in zip(ra, dec, psi, rnd):
		if n < sky_params(net, time, r, d, p):
		ndraws -= 1
		ra_out[ndraws] = r
		dec_out[ndraws] = d
		psi_out[ndraws] = p
		break
		return ra_out, dec_out, psi_out

		def parse_polarisation(self, polarisation):
		"""
		@@ -117,38 +81,72 @@ class Waveform(object):
		plt.plot(self.hp.data.data, label="+ polarisation")
		plt.plot(self.hx.data.data, label="x polarisation")


		def uniform_time(self):
		def _generate(self, rate=16384.0):
		"""
		Get a set of randomized (integer) event times.
		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
		"""
		return random.randint(self.tstart, self.tstop, self.expnum) + random.rand(self.expnum)
		row = self._row()
		self.swig_row = lalburst.CreateSimBurst()
		for a in lsctables.SimBurstTable.validcolumns.keys():
		try:
		setattr(self.swig_row, a, getattr( row, a ))
		print getattr(row,a)
		except AttributeError: continue
		except TypeError:
		continue
		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):
		def _row(self, sim=None):
		"""
		Produce a simburst table row for this waveform.

		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.
		"""
		if not sim: sim = self.sim
		self.row = sim.RowType()

		# Required columns not defined makes ligolw unhappy
		for a in lsctables.SimBurstTable.validcolumns.keys():
		setattr(self.row, a, None)
		setattr(self.row, a, self.params[a])

		self.row.waveform = self.waveform
		# Fill in the 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()
		# Get the sky locations
		self.row.ra, self.row.dec, self.row.psi = self.sky_dist()
		self.row.simulation_id = sim.get_next_id()
		self.row.waveform_number = random.randint(0,int(2**32)-1)
		### !! 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)
		self.row.time_slide_id = ilwd.ilwdchar("time_slide:time_slide_id:%d" % 1)#options.time_slide_id)

		return self.row



		@@ -158,178 +156,90 @@ class SineGaussian(Waveform):
		"""
		waveform = "SineGaussian"

		def __init__(self, q, frequency, hrss, polarisation, time, seed=0):
		def __init__(self, q, frequency, hrss, polarisation, time, sky_dist=uniform_sky, seed=0):
		"""A class to represent a SineGaussian ad-hoc waveform.

		Parameters
		----------
		q : float or list
		q : float
		The quality factor.
		If a float is provded then the q-factor is fixed. If a list is provided then they are the
		minimum and maximum quality factor of the injections.

		frequency : float or list
		frequency : float
		The frequency of the injection.
		If a float is provided the frequency will be fixed, if a list is provided then this will be the
		minimum and maximum frequencies.

		hrss : float or list
		hrss : float
		The strain magnitude of the injection.
		If a float is provided then the hrss will be fixed, if a list is provided then this will be the
		minimum and maximum hrss.

		polarisation : str {'linear', 'elliptical', 'circular'}
		The type of polarisation of the waveform.

		time : float or list
		The time period over which the injection should be made. If
		a list is given they should be the start and end times, and
		the waveform will be produced at some random point in that
		time range. If a float is given then the injection will be
		made at that specific time.
		time : float
		The central time of the injection.

		seed : float
		sky_dist : func
		The function describing the sky distribution which the injections
		should be made over. Defaults to a uniform sky.

		seed : int
		The random seed used to make the injection time of the waveform.
		The default seed is 0.

		"""

		if isinstance(q, list):
		self.q_values = segment(q[0], q[1])
		else:
		self.q_values = q

		if isinstance(frequency, list):
		self.f_values = segment(frequency[0], frequency[1])
		else:
		self.f_values = frequency

		if isinstance(hrss, list):
		self.h_values = segment(hrss[0], hrss[1])
		else:
		self.h_values = hrss

		if isinstance(time, list):
		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))[0]
		else:
		self._clear_params()
		self.sky_dist = sky_dist
		self.params['hrss'] = hrss
		self.params['frequency'] = frequency
		self.params['q'] = q
		self.time = time

		self.polarisation = polarisation
		self.params['pol_ellipse_e'], self.params['pol_ellipse_angle'] = self.parse_polarisation(self.polarisation)

		self.q, self.frequency, self.hrss = self.draw_params()

		def draw_params(self):
		class Gaussian(Waveform):
		"""
		Draw a set of parameters (hrss, q, frequency) appropriate for a set of sinegaussian type injections.

		Returns
		-------
		q : float
		Quality factor
		f : float
		Frequency
		h : float
		HRSS
		A class to represent a Gaussian injection.
		"""
		h, q, f = [], [], []
		if type(self.q_values) == segment:
		q = uniform_interval(self.q_valies, 1)
		else:
		q = self.q_values
		#q = map(lambda i: self.q_range[i], random.randint(0, len(self.q_range), self.expnum) )

		if type(self.h_values) == segment:
		h = self.log_hrss()
		else:
		h = self.h_values
		#h = map(lambda i: self.h_range[i], random.randint(0, len(self.h_range), self.expnum) )

		if type(self.f_values) == segment:
		f = uniform_interval(self.f_values, self.expnum)
		else:
		f = self.f_values
		#f = map(lambda i: self.f_range[i], random.randint(0, len(self.f_range), self.expnum))

		return q, f, h

		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
		waveform = "Gaussian"

		def _row(self, sim=None):
		def __init__(self, duration, hrss, time, sky_dist=uniform_sky, seed=0):
		"""
		Produce a simburst table row for this waveform.
		A class to represent a Gaussian ad-hoc waveform.

		Parameters
		----------
		sim : lsctable
		The sim table for this row.
		duration : float or list
		The duration, in milliseconds, of the Gaussian waveform.

		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(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 = "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)
		hrss : float or list
		The strain magnitude of the injection.
		If a float is provided then the hrss will be fixed,
		if a list is provided then this will be the
		minimum and maximum hrss.

		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)
		polarisation : str {'linear', 'elliptical', 'circular'}
		The type of polarisation of the waveform.

		return self.row
		time : float or list
		The time period over which the injection should be made. If
		a list is given they should be the start and end times, and
		the waveform will be produced at some random point in that
		time range. If a float is given then the injection will be
		made at that specific time.

		#def _repr_html_(self):
		sky_dist : func
		The function describing the sky distribution which the injections
		should be made over. Defaults to a uniform sky.

		seed : float
		The random seed used to make the injection time of the waveform.
		The default seed is 0.

		"""
		self._clear_params()
		self.sky_dist = sky_dist
		self.params['duration'] = duration
		self.params['hrss'] = hrss
		self.time = time
		self.params['pol_ellipse_e'] = 1.0
		self.params['pol_ellipse_angle'] = 0

notebooks/Untitled.ipynb

deleted100644 → 0

+0 −241

File deleted.

Preview size limit exceeded, changes collapsed.