inspiral_lr.py 43.1 KB
Newer Older
Kipp Cannon's avatar
Kipp Cannon committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
# Copyright (C) 2017  Kipp Cannon
# Copyright (C) 2011--2014  Kipp Cannon, Chad Hanna, Drew Keppel
# Copyright (C) 2013  Jacob Peoples
#
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 2 of the License, or (at your
# option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.


#
# =============================================================================
#
#                                   Preamble
#
# =============================================================================
#


29 30 31 32 33
try:
	from fpconst import NegInf
except ImportError:
	# not all machines have fpconst installed
	NegInf = float("-inf")
Kipp Cannon's avatar
Kipp Cannon committed
34 35
import math
import numpy
36
import os
Kipp Cannon's avatar
Kipp Cannon committed
37 38 39 40
import random
from scipy import interpolate
from scipy import stats
import sys
41
import warnings
42
import json
Kipp Cannon's avatar
Kipp Cannon committed
43 44


45 46 47 48 49
from ligo.lw import ligolw
from ligo.lw import lsctables
from ligo.lw import array as ligolw_array
from ligo.lw import param as ligolw_param
from ligo.lw import utils as ligolw_utils
50
from ligo import segments
Kipp Cannon's avatar
Kipp Cannon committed
51 52
from gstlal.stats import horizonhistory
from gstlal.stats import inspiral_extrinsics
53
from gstlal.stats import inspiral_intrinsics
Kipp Cannon's avatar
Kipp Cannon committed
54 55 56 57 58 59 60
from gstlal.stats import trigger_rate
import lal
from lal import rate
from lalburst import snglcoinc
import lalsimulation


61 62 63 64 65 66 67
# FIXME:  caution, this information might get organized differently later.
# for now we just need to figure out where the gstlal-inspiral directory in
# share/ is.  don't write anything that assumes that this module will
# continue to define any of these symbols
from gstlal import paths as gstlal_config_paths


Kipp Cannon's avatar
Kipp Cannon committed
68 69
__all__ = [
	"LnSignalDensity",
70 71 72 73 74
	"LnNoiseDensity",
	"DatalessLnSignalDensity",
	"DatalessLnNoiseDensity",
	"OnlineFrankensteinLnSignalDensity",
	"OnlineFrankensteinLnNoiseDensity"
Kipp Cannon's avatar
Kipp Cannon committed
75 76 77
]


78 79 80 81 82
# The typical horizon distance used to normalize such values in the likelihood
# ratio. Units are Mpc
TYPICAL_HORIZON_DISTANCE = 150.


Kipp Cannon's avatar
Kipp Cannon committed
83 84 85 86 87 88 89 90 91 92 93
#
# =============================================================================
#
#                                  Base Class
#
# =============================================================================
#


class LnLRDensity(snglcoinc.LnLRDensity):
	# range of SNRs covered by this object
94
	snr_min = 4.0
Kipp Cannon's avatar
Kipp Cannon committed
95 96 97 98

	# SNR, \chi^2 binning definition
	snr_chi_binning = rate.NDBins((rate.ATanLogarithmicBins(2.6, 26., 300), rate.ATanLogarithmicBins(.001, 0.2, 280)))

99
	def __init__(self, template_ids, instruments, delta_t, min_instruments = 2):
Kipp Cannon's avatar
Kipp Cannon committed
100 101 102 103
		#
		# check input
		#

104 105
		if template_ids is not None and not template_ids:
			raise ValueError("template_ids cannot be empty")
Kipp Cannon's avatar
Kipp Cannon committed
106 107 108 109 110 111 112 113 114 115 116
		if min_instruments < 1:
			raise ValueError("min_instruments=%d must be >= 1" % min_instruments)
		if min_instruments > len(instruments):
			raise ValueError("not enough instruments (%s) to satisfy min_instruments=%d" % (", ".join(sorted(instruments)), min_instruments))
		if delta_t < 0.:
			raise ValueError("delta_t=%g must be >= 0" % delta_t)

		#
		# initialize
		#

117
		self.template_ids = frozenset(template_ids) if template_ids is not None else template_ids
Kipp Cannon's avatar
Kipp Cannon committed
118 119 120 121 122 123 124 125 126 127 128
		self.instruments = frozenset(instruments)
		self.min_instruments = min_instruments
		self.delta_t = delta_t
		# install SNR, chi^2 PDFs.  numerator will override this
		self.densities = {}
		for instrument in instruments:
			self.densities["%s_snr_chi" % instrument] = rate.BinnedLnPDF(self.snr_chi_binning)

	def __iadd__(self, other):
		if type(other) != type(self):
			raise TypeError(other)
129 130 131 132 133 134
		# template_id set mismatch is allowed in the special case
		# that one or the other is None to make it possible to
		# construct generic seed objects providing initialization
		# data for the ranking statistics.
		if self.template_ids is not None and other.template_ids is not None and self.template_ids != other.template_ids:
			raise ValueError("incompatible template IDs")
Kipp Cannon's avatar
Kipp Cannon committed
135 136 137 138 139 140 141
		if self.instruments != other.instruments:
			raise ValueError("incompatible instrument sets")
		if self.min_instruments != other.min_instruments:
			raise ValueError("incompatible minimum number of instruments")
		if self.delta_t != other.delta_t:
			raise ValueError("incompatible delta_t coincidence thresholds")

142 143 144
		if self.template_ids is None and other.template_ids is not None:
			self.template_ids = frozenset(other.template_ids)

Kipp Cannon's avatar
Kipp Cannon committed
145 146 147 148 149 150 151 152 153 154
		for key, lnpdf in self.densities.items():
			lnpdf += other.densities[key]

		try:
			del self.interps
		except AttributeError:
			pass
		return self

	def increment(self, event):
155 156 157
		# this test is intended to fail if .template_ids is None:
		# must not collect trigger statistics unless we can verify
		# that they are for the correct templates.
158
		if event.template_id not in self.template_ids:
159
			raise ValueError("event from wrong template")
160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
		# ignore events below threshold.  the trigger generation
		# mechanism now produces sub-threshold triggers to
		# facilitiate the construction of sky maps in the online
		# search.  it's therefore now our responsibility to exclude
		# them from the PDFs in the ranking statistic here.
		#
		# FIXME:  this creates an inconsistency.  we exclude events
		# from the PDFs, but we don't exclude them from the
		# evaluation of those PDFs in the .__call__() method,
		# instead we'll report that the candidate is impossible
		# because it contains a trigger in a "forbidden" region of
		# the parameter space (forbidden because we've ensured
		# nothing populates those bins).  the .__call__() methods
		# should probably be modified to exclude triggers with
		# sub-threshold SNRs just as we've excluded them here, but
		# I'm not yet sure.  instead, for the time being, we
		# resolve the problem by excluding sub-threshold triggers
		# from the "-from-triggers" method in far.py.  there is a
		# FIXME related to this same issue in that code.
		if event.snr < self.snr_min:
			return
Kipp Cannon's avatar
Kipp Cannon committed
181 182 183
		self.densities["%s_snr_chi" % event.ifo].count[event.snr, event.chisq / event.snr**2.] += 1.0

	def copy(self):
184
		new = type(self)(self.template_ids, self.instruments, self.delta_t, self.min_instruments)
Kipp Cannon's avatar
Kipp Cannon committed
185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247
		for key, lnpdf in self.densities.items():
			new.densities[key] = lnpdf.copy()
		return new

	def mkinterps(self):
		self.interps = dict((key, lnpdf.mkinterp()) for key, lnpdf in self.densities.items())

	def finish(self):
		snr_kernel_width_at_8 = 8.,
		chisq_kernel_width = 0.08,
		sigma = 10.
		for key, lnpdf in self.densities.items():
			# construct the density estimation kernel.  be
			# extremely conservative and assume only 1 in 10
			# samples are independent, but assume there are
			# always at least 1e7 samples.
			numsamples = max(lnpdf.array.sum() / 10. + 1., 1e6)
			snr_bins, chisq_bins = lnpdf.bins
			snr_per_bin_at_8 = (snr_bins.upper() - snr_bins.lower())[snr_bins[8.]]
			chisq_per_bin_at_0_02 = (chisq_bins.upper() - chisq_bins.lower())[chisq_bins[0.02]]

			# apply Silverman's rule so that the width scales
			# with numsamples**(-1./6.) for a 2D PDF
			snr_kernel_bins = snr_kernel_width_at_8 / snr_per_bin_at_8 / numsamples**(1./6.)
			chisq_kernel_bins = chisq_kernel_width / chisq_per_bin_at_0_02 / numsamples**(1./6.)

			# check the size of the kernel. We don't ever let
			# it get smaller than the 2.5 times the bin size
			if snr_kernel_bins < 2.5:
				snr_kernel_bins = 2.5
				warnings.warn("Replacing snr kernel bins with 2.5")
			if chisq_kernel_bins < 2.5:
				chisq_kernel_bins = 2.5
				warnings.warn("Replacing chisq kernel bins with 2.5")

			# convolve bin count with density estimation kernel
			rate.filter_array(lnpdf.array, rate.gaussian_window(snr_kernel_bins, chisq_kernel_bins, sigma = sigma))

			# zero everything below the SNR cut-off.  need to
			# do the slicing ourselves to avoid zeroing the
			# at-threshold bin
			lnpdf.array[:lnpdf.bins[0][self.snr_min],:] = 0.

			# normalize what remains
			lnpdf.normalize()
		self.mkinterps()

		#
		# never allow PDFs that have had the density estimation
		# transform applied to be written to disk:  on-disk files
		# must only ever provide raw counts.  also don't allow
		# density estimation to be applied twice
		#

		def to_xml(*args, **kwargs):
			raise NotImplementedError("writing .finish()'ed LnLRDensity object to disk is forbidden")
		self.to_xml = to_xml
		def finish(*args, **kwargs):
			raise NotImplementedError(".finish()ing a .finish()ed LnLRDensity object is forbidden")
		self.finish = finish

	def to_xml(self, name):
		xml = super(LnLRDensity, self).to_xml(name)
248
		xml.appendChild(ligolw_param.Param.from_pyvalue(u"template_ids", ",".join("%d" % template_id for template_id in sorted(self.template_ids)) if self.template_ids is not None else None))
249
		# FIXME this is not an ideal way to get only one into the file
Kipp Cannon's avatar
Kipp Cannon committed
250 251 252 253 254 255 256 257 258 259 260 261 262
		xml.appendChild(ligolw_param.Param.from_pyvalue(u"instruments", lsctables.instrumentsproperty.set(self.instruments)))
		xml.appendChild(ligolw_param.Param.from_pyvalue(u"min_instruments", self.min_instruments))
		xml.appendChild(ligolw_param.Param.from_pyvalue(u"delta_t", self.delta_t))
		for key, lnpdf in self.densities.items():
			# never write PDFs to disk without ensuring the
			# normalization metadata is up to date
			lnpdf.normalize()
			xml.appendChild(lnpdf.to_xml(key))
		return xml

	@classmethod
	def from_xml(cls, xml, name):
		xml = cls.get_xml_root(xml, name)
263
		template_ids = ligolw_param.get_pyvalue(xml, u"template_ids")
264 265 266 267 268 269 270
		# FIXME find a better way to do this: Allow a file to not have
		# a mass model filename.  This is can happen when e.g., gstlal
		# inspiral produces a ranking stat file which just records but
		# doesn't evaluate LRs.  Then you marginalize it in with a
		# create prior diststats file which does have the mass model.
		# The iadd method makes sure that the mass model is unique or
		# that it doesn't exist.
271 272
		if template_ids is not None:
			template_ids = frozenset(int(i) for i in template_ids.split(","))
Kipp Cannon's avatar
Kipp Cannon committed
273
		self = cls(
274
			template_ids = template_ids,
Kipp Cannon's avatar
Kipp Cannon committed
275 276
			instruments = lsctables.instrumentsproperty.get(ligolw_param.get_pyvalue(xml, u"instruments")),
			delta_t = ligolw_param.get_pyvalue(xml, u"delta_t"),
277
			min_instruments = ligolw_param.get_pyvalue(xml, u"min_instruments")
Kipp Cannon's avatar
Kipp Cannon committed
278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293
		)
		for key in self.densities:
			self.densities[key] = self.densities[key].from_xml(xml, key)
		return self


#
# =============================================================================
#
#                              Numerator Density
#
# =============================================================================
#


class LnSignalDensity(LnLRDensity):
294
	def __init__(self, *args, **kwargs):
Kipp Cannon's avatar
Kipp Cannon committed
295
		population_model_file = kwargs.pop("population_model_file", None)
296
		dtdphi_file = kwargs.pop("dtdphi_file", None)
297
		self.horizon_factors = kwargs.pop("horizon_factors", None)
298
		super(LnSignalDensity, self).__init__(*args, **kwargs)
Kipp Cannon's avatar
Kipp Cannon committed
299 300 301 302 303 304 305
		# install SNR, chi^2 PDF (one for all instruments)
		self.densities = {
			"snr_chi": inspiral_extrinsics.NumeratorSNRCHIPDF(self.snr_chi_binning)
		}

		# record of horizon distances for all instruments in the
		# network
306
		self.horizon_history = horizonhistory.HorizonHistories((instrument, horizonhistory.NearestLeafTree()) for instrument in self.instruments)
Kipp Cannon's avatar
Kipp Cannon committed
307
		self.population_model_file = population_model_file
308
		self.dtdphi_file = dtdphi_file
309 310
		# source population model
		# FIXME:  introduce a mechanism for selecting the file
311
		self.population_model = inspiral_intrinsics.SourcePopulationModel(self.template_ids, filename = self.population_model_file)
312 313
		if self.dtdphi_file is not None:
			self.InspiralExtrinsics = inspiral_extrinsics.InspiralExtrinsics(self.min_instruments, self.dtdphi_file)
314

315 316 317
	def set_horizon_factors(self, horizon_factors):
		self.horizon_factors = horizon_factors

318
	def __call__(self, segments, snrs, chi2s_over_snr2s, phase, dt, template_id):
Kipp Cannon's avatar
Kipp Cannon committed
319 320
		assert frozenset(segments) == self.instruments
		if len(snrs) < self.min_instruments:
321
			return NegInf
Kipp Cannon's avatar
Kipp Cannon committed
322

323 324
		# use volume-weighted average horizon distance over
		# duration of event to estimate sensitivity
Kipp Cannon's avatar
Kipp Cannon committed
325 326 327
		assert all(segments.values()), "encountered trigger with duration = 0"
		horizons = dict((instrument, (self.horizon_history[instrument].functional_integral(map(float, seg), lambda d: d**3.) / float(abs(seg)))**(1./3.)) for instrument, seg in segments.items())

328 329 330 331 332 333
		# compute P(t).  P(t) \propto volume within which a signal will
		# produce a candidate * number of trials \propto cube of
		# distance to which the mininum required number of instruments
		# can see (I think) * number of templates.  we measure the
		# distance in multiples of TYPICAL_HORIZON_DISTANCE Mpc just to
		# get a number that will be, typically, a little closer to 1,
Kipp Cannon's avatar
Kipp Cannon committed
334
		# in lieu of properly normalizating this factor.  we can't
335 336 337 338 339 340 341 342 343 344 345 346 347
		# normalize the factor because the normalization depends on the
		# duration of the experiment, and that keeps changing when
		# running online, combining offline chunks from different
		# times, etc., and that would prevent us from being able to
		# compare a ln L ranking statistic value defined during one
		# period to ranking statistic values defined in other periods.
		# by removing the normalization, and leaving this be simply a
		# factor that is proportional to the rate of signals, we can
		# compare ranking statistics across analysis boundaries but we
		# loose meaning in the overall scale:  ln L = 0 is not a
		# special value, as it would be if the numerator and
		# denominator were both normalized properly.
		horizon = sorted(horizons.values())[-self.min_instruments] / TYPICAL_HORIZON_DISTANCE
348
		if not horizon:
349
			return NegInf
350 351 352 353 354 355
		# horizon factors adjusts the computed horizon factor according
		# to the sigma values computed at the time of the SVD. This
		# gives a good approximation to the horizon distance for each
		# template without computing them each explicitly. Only one
		# template has its horizon calculated explicitly.
		lnP = 3. * math.log(horizon * self.horizon_factors[template_id]) + math.log(len(self.template_ids))
356

357 358 359 360 361
		# Add P(instruments | horizon distances)
		try:
			lnP += math.log(self.InspiralExtrinsics.p_of_instruments_given_horizons(snrs.keys(), horizons))
		except ValueError:
			# The code raises a value error when a needed horizon distance is zero
362
			return NegInf
Kipp Cannon's avatar
Kipp Cannon committed
363

364 365 366 367 368
		# Evaluate dt, dphi, snr probability
		try:
			lnP += math.log(self.InspiralExtrinsics.time_phase_snr(dt, phase, snrs, horizons))
		# FIXME need to make sure this is really a math domain error
		except ValueError:
369
			return NegInf
Kipp Cannon's avatar
Kipp Cannon committed
370

371 372 373
		# evaluate population model
		lnP += self.population_model.lnP_template_signal(template_id, max(snrs.values()))

374 375 376
		# evalute the (snr, \chi^2 | snr) PDFs (same for all
		# instruments)
		interp = self.interps["snr_chi"]
377
		return lnP + sum(interp(snrs[instrument], chi2_over_snr2) for instrument, chi2_over_snr2 in chi2s_over_snr2s.items())
Kipp Cannon's avatar
Kipp Cannon committed
378 379 380 381

	def __iadd__(self, other):
		super(LnSignalDensity, self).__iadd__(other)
		self.horizon_history += other.horizon_history
382 383 384 385
		if self.population_model_file is not None and other.population_model_file is not None and other.population_model_file != self.population_model_file:
			raise ValueError("incompatible mass model file names")
		if self.population_model_file is None and other.population_model_file is not None:
			self.population_model_file = other.population_model_file
386 387 388 389
		if self.dtdphi_file is not None and other.dtdphi_file is not None and other.dtdphi_file != self.dtdphi_file:
			raise ValueError("incompatible dtdphi files")
		if self.dtdphi_file is None and other.dtdphi_file is not None:
			self.dtdphi_file = other.dtdphi_file
390
		if self.horizon_factors is not None and other.horizon_factors is not None and other.horizon_factors != self.horizon_factors:
391 392 393 394
			# require that the horizon factors be the same within 1%
			for k in self.horizon_factors:
				if 0.99 > self.horizon_factors[k] / other.horizon_factors[k] > 1.01:
					raise ValueError("incompatible horizon_factors")
395 396
		if self.horizon_factors is None and other.horizon_factors is not None:
			self.horizon_factors = other.horizon_factors
397

Kipp Cannon's avatar
Kipp Cannon committed
398 399 400 401 402 403 404 405
		return self

	def increment(self, *args, **kwargs):
		raise NotImplementedError

	def copy(self):
		new = super(LnSignalDensity, self).copy()
		new.horizon_history = self.horizon_history.copy()
Kipp Cannon's avatar
Kipp Cannon committed
406
		new.population_model_file = self.population_model_file
407
		new.dtdphi_file = self.dtdphi_file
Kipp Cannon's avatar
Kipp Cannon committed
408 409 410
		# okay to use references because read-only data
		new.population_model = self.population_model
		new.InspiralExtrinsics = self.InspiralExtrinsics
411
		new.horizon_factors = self.horizon_factors
Kipp Cannon's avatar
Kipp Cannon committed
412 413 414 415 416 417 418 419 420 421 422 423 424 425
		return new

	def local_mean_horizon_distance(self, gps, window = segments.segment(-32., +2.)):
		# horizon distance window is in seconds.  this is sort of a
		# hack, we should really tie this to each waveform's filter
		# length somehow, but we don't have a way to do that and
		# it's not clear how much would be gained by going to the
		# trouble of implementing that.  I don't even know what to
		# set it to, so for now it's set to something like the
		# typical width of the whitening filter's impulse response.
		t = abs(window)
		vtdict = self.horizon_history.functional_integral_dict(window.shift(float(gps)), lambda D: D**3.)
		return dict((instrument, (vt / t)**(1./3.)) for instrument, vt in vtdict.items())

426
	def add_signal_model(self, prefactors_range = (0.001, 0.30), df = 200, inv_snr_pow = 4.):
Kipp Cannon's avatar
Kipp Cannon committed
427 428
		# normalize to 10 *mi*llion signals.  this count makes the
		# density estimation code choose a suitable kernel size
429
		inspiral_extrinsics.NumeratorSNRCHIPDF.add_signal_model(self.densities["snr_chi"], 1e12, prefactors_range, df, inv_snr_pow = inv_snr_pow, snr_min = self.snr_min)
Kipp Cannon's avatar
Kipp Cannon committed
430 431 432 433 434 435 436 437 438 439 440 441 442 443
		self.densities["snr_chi"].normalize()

	def candidate_count_model(self, rate = 1000.):
		"""
		Compute and return a prediction for the total number of
		above-threshold signal candidates expected.  The rate
		parameter sets the nominal signal rate in units of Gpc^-3
		a^-1.
		"""
		# FIXME:  this needs to understand a mass distribution
		# model and what part of the mass space this numerator PDF
		# is for
		seg = (self.horizon_history.minkey(), self.horizon_history.maxkey())
		V_times_t = self.horizon_history.functional_integral(seg, lambda horizons: sorted(horizons.values())[-self.min_instruments]**3.)
444 445 446
		# Mpc**3 --> Gpc**3, seconds --> years
		V_times_t *= 1e-9 / (86400. * 365.25)
		return V_times_t * rate * len(self.template_ids)
Kipp Cannon's avatar
Kipp Cannon committed
447 448 449 450 451 452 453

	def random_sim_params(self, sim, horizon_distance = None, snr_efficiency = 1.0, coinc_only = True):
		"""
		Generator that yields an endless sequence of randomly
		generated parameter dictionaries drawn from the
		distribution of parameters expected for the given
		injection, which is an instance of a SimInspiral table row
454
		object (see ligo.lw.lsctables.SimInspiral for more
Kipp Cannon's avatar
Kipp Cannon committed
455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480
		information).  Each value in the sequence is a tuple, the
		first element of which is the random parameter dictionary
		and the second is 0.

		See also:

		LnNoiseDensity.random_params()

		The sequence is suitable for input to the .ln_lr_samples()
		log likelihood ratio generator.

		Bugs:

		The second element in each tuple in the sequence is merely
		a placeholder, not the natural logarithm of the PDF from
		which the sample has been drawn, as in the case of
		random_params().  Therefore, when used in combination with
		.ln_lr_samples(), the two probability densities computed
		and returned by that generator along with each log
		likelihood ratio value will simply be the probability
		densities of the signal and noise populations at that point
		in parameter space.  They cannot be used to form an
		importance weighted sampler of the log likelihood ratios.
		"""
		# FIXME:  this is still busted since the rewrite

481
		# FIXME need to add dt and dphi
Kipp Cannon's avatar
Kipp Cannon committed
482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553
		#
		# retrieve horizon distance from history if not given
		# explicitly.  retrieve SNR threshold from class attribute
		# if not given explicitly
		#

		if horizon_distance is None:
			horizon_distance = self.local_mean_horizon_distance(sim.time_geocent)

		#
		# compute nominal SNRs
		#

		cosi2 = math.cos(sim.inclination)**2.
		gmst = lal.GreenwichMeanSiderealTime(sim.time_geocent)
		snr_0 = {}
		for instrument, DH in horizon_distance.items():
			fp, fc = lal.ComputeDetAMResponse(lalsimulation.DetectorPrefixToLALDetector(str(instrument)).response, sim.longitude, sim.latitude, sim.polarization, gmst)
			snr_0[instrument] = snr_efficiency * 8. * DH * math.sqrt(fp**2. * (1. + cosi2)**2. / 4. + fc**2. * cosi2) / sim.distance

		#
		# construct SNR generators, and approximating the SNRs to
		# be fixed at the nominal SNRs construct \chi^2 generators
		#

		def snr_gen(snr):
			rvs = stats.ncx2(2., snr**2.).rvs
			math_sqrt = math.sqrt
			while 1:
				yield math_sqrt(rvs())

		def chi2_over_snr2_gen(instrument, snr):
			rates_lnx = numpy.log(self.injection_rates["%s_snr_chi" % instrument].bins[1].centres())
			# FIXME:  kinda broken for SNRs below self.snr_min
			rates_cdf = self.injection_rates["%s_snr_chi" % instrument][max(snr, self.snr_min),:].cumsum()
			# add a small tilt to break degeneracies then
			# normalize
			rates_cdf += numpy.linspace(0., 0.001 * rates_cdf[-1], len(rates_cdf))
			rates_cdf /= rates_cdf[-1]
			assert not numpy.isnan(rates_cdf).any()

			interp = interpolate.interp1d(rates_cdf, rates_lnx)
			math_exp = math.exp
			random_random = random.random
			while 1:
				yield math_exp(float(interp(random_random())))

		gens = dict(((instrument, "%s_snr_chi" % instrument), (iter(snr_gen(snr)).next, iter(chi2_over_snr2_gen(instrument, snr)).next)) for instrument, snr in snr_0.items())

		#
		# yield a sequence of randomly generated parameters for
		# this sim.
		#

		while 1:
			params = {"snrs": {}}
			instruments = []
			for (instrument, key), (snr, chi2_over_snr2) in gens.items():
				snr = snr()
				if snr < self.snr_min:
					continue
				params[key] = snr, chi2_over_snr2()
				params["snrs"][instrument] = snr
				instruments.append(instrument)
			if coinc_only and len(instruments) < self.denominator.min_instruments:
				continue
			params.horizons = horizon_distance
			yield params, 0.

	def to_xml(self, name):
		xml = super(LnSignalDensity, self).to_xml(name)
		xml.appendChild(self.horizon_history.to_xml(u"horizon_history"))
554
		xml.appendChild(ligolw_param.Param.from_pyvalue(u"population_model_file", self.population_model_file))
555
		xml.appendChild(ligolw_param.Param.from_pyvalue(u"horizon_factors", json.dumps(self.horizon_factors) if self.horizon_factors is not None else None))
556
		xml.appendChild(ligolw_param.Param.from_pyvalue(u"dtdphi_file", self.dtdphi_file))
Kipp Cannon's avatar
Kipp Cannon committed
557 558 559 560 561 562 563
		return xml

	@classmethod
	def from_xml(cls, xml, name):
		xml = cls.get_xml_root(xml, name)
		self = super(LnSignalDensity, cls).from_xml(xml, name)
		self.horizon_history = horizonhistory.HorizonHistories.from_xml(xml, u"horizon_history")
Kipp Cannon's avatar
Kipp Cannon committed
564
		self.population_model_file = ligolw_param.get_pyvalue(xml, u"population_model_file")
565
		self.dtdphi_file = ligolw_param.get_pyvalue(xml, u"dtdphi_file")
566 567 568 569 570 571 572
		self.horizon_factors = ligolw_param.get_pyvalue(xml, u"horizon_factors")
		if self.horizon_factors is not None:
			# FIXME, how do we properly decode the json, I assume something in ligolw can do this?
			self.horizon_factors = self.horizon_factors.replace("\\","").replace('\\"','"')
			self.horizon_factors = json.loads(self.horizon_factors)
			self.horizon_factors = dict((int(k), v) for k, v in self.horizon_factors.items())
			assert set(self.template_ids) == set(self.horizon_factors)
Kipp Cannon's avatar
Kipp Cannon committed
573
		self.population_model = inspiral_intrinsics.SourcePopulationModel(self.template_ids, filename = self.population_model_file)
574 575
		if self.dtdphi_file is not None:
			self.InspiralExtrinsics = inspiral_extrinsics.InspiralExtrinsics(self.min_instruments, self.dtdphi_file)
Kipp Cannon's avatar
Kipp Cannon committed
576 577 578 579 580 581 582 583 584 585 586
		return self


class DatalessLnSignalDensity(LnSignalDensity):
	"""
	Stripped-down version of LnSignalDensity for use in estimating
	ranking statistics when no data has been collected from an
	instrument with which to define the ranking statistic.

	Used, for example, to implement low-significance candidate culls,
	etc.
587 588 589 590

	Assumes all available instruments are on and have the same horizon
	distance, and assess candidates based only on SNR and \chi^2
	distributions.
Kipp Cannon's avatar
Kipp Cannon committed
591
	"""
592 593 594 595
	def __init__(self, *args, **kwargs):
		super(DatalessLnSignalDensity, self).__init__(*args, **kwargs)
		# so we're ready to go!
		self.add_signal_model()
Kipp Cannon's avatar
Kipp Cannon committed
596

597
	def __call__(self, segments, snrs, chi2s_over_snr2s, phase, dt, template_id):
598 599
		# evaluate P(t) \propto number of templates
		lnP = math.log(len(self.template_ids))
Kipp Cannon's avatar
Kipp Cannon committed
600

Kipp Cannon's avatar
Kipp Cannon committed
601 602 603 604
		# Add P(instruments | horizon distances).  assume all
		# instruments have TYPICAL_HORIZON_DISTANCE horizon
		# distance
		horizons = dict.fromkeys(segments, TYPICAL_HORIZON_DISTANCE)
605 606 607
		try:
			lnP += math.log(self.InspiralExtrinsics.p_of_instruments_given_horizons(snrs.keys(), horizons))
		except ValueError:
Kipp Cannon's avatar
Kipp Cannon committed
608 609
			# raises ValueError when a needed horizon distance
			# is zero
610
			return NegInf
611 612 613 614 615 616

		# Evaluate dt, dphi, snr probability
		try:
			lnP += math.log(self.InspiralExtrinsics.time_phase_snr(dt, phase, snrs, horizons))
		# FIXME need to make sure this is really a math domain error
		except ValueError:
617
			return NegInf
Kipp Cannon's avatar
Kipp Cannon committed
618

619 620 621
		# evaluate population model
		lnP += self.population_model.lnP_template_signal(template_id, max(snrs.values()))

622 623 624
		# evalute the (snr, \chi^2 | snr) PDFs (same for all
		# instruments)
		interp = self.interps["snr_chi"]
625
		return lnP + sum(interp(snrs[instrument], chi2_over_snr2) for instrument, chi2_over_snr2 in chi2s_over_snr2s.items())
Kipp Cannon's avatar
Kipp Cannon committed
626 627 628 629

	def __iadd__(self, other):
		raise NotImplementedError

630 631 632
	def increment(self, *args, **kwargs):
		raise NotImplementedError

Kipp Cannon's avatar
Kipp Cannon committed
633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649
	def copy(self):
		raise NotImplementedError

	def to_xml(self, name):
		# I/O not permitted:  the on-disk version would be
		# indistinguishable from a real ranking statistic and could
		# lead to accidents
		raise NotImplementedError

	@classmethod
	def from_xml(cls, xml, name):
		# I/O not permitted:  the on-disk version would be
		# indistinguishable from a real ranking statistic and could
		# lead to accidents
		raise NotImplementedError


Kipp Cannon's avatar
Kipp Cannon committed
650
class OnlineFrankensteinLnSignalDensity(LnSignalDensity):
651 652 653 654 655 656 657 658 659 660 661 662 663 664
	"""
	Version of LnSignalDensity with horizon distance history spliced in
	from another instance.  Used to solve a chicken-or-egg problem and
	assign ranking statistic values in an aonline anlysis.  NOTE:  the
	horizon history is not copied from the donor, instances of this
	class hold a reference to the donor's data, so as it is modified
	those modifications are immediately reflected here.

	For safety's sake, instances cannot be written to or read from
	files, cannot be marginalized together with other instances, nor
	accept updates from new data.
	"""
	@classmethod
	def splice(cls, src, Dh_donor):
665
		self = cls(src.template_ids, src.instruments, src.delta_t, population_model_file = src.population_model_file, min_instruments = src.min_instruments, horizon_factors = src.horizon_factors)
666
		for key, lnpdf in src.densities.items():
667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695
			self.densities[key] = lnpdf.copy()
		# NOTE:  not a copy.  we hold a reference to the donor's
		# data so that as it is updated, we get the updates.
		self.horizon_history = Dh_donor.horizon_history
		return self

	def __iadd__(self, other):
		raise NotImplementedError

	def increment(self, *args, **kwargs):
		raise NotImplementedError

	def copy(self):
		raise NotImplementedError

	def to_xml(self, name):
		# I/O not permitted:  the on-disk version would be
		# indistinguishable from a real ranking statistic and could
		# lead to accidents
		raise NotImplementedError

	@classmethod
	def from_xml(cls, xml, name):
		# I/O not permitted:  the on-disk version would be
		# indistinguishable from a real ranking statistic and could
		# lead to accidents
		raise NotImplementedError


Kipp Cannon's avatar
Kipp Cannon committed
696 697 698 699 700 701 702 703 704 705
#
# =============================================================================
#
#                             Denominator Density
#
# =============================================================================
#


class LnNoiseDensity(LnLRDensity):
706 707
	def __init__(self, *args, **kwargs):
		super(LnNoiseDensity, self).__init__(*args, **kwargs)
Kipp Cannon's avatar
Kipp Cannon committed
708 709 710

		# record of trigger counts vs time for all instruments in
		# the network
711
		self.triggerrates = trigger_rate.triggerrates((instrument, trigger_rate.ratebinlist()) for instrument in self.instruments)
Kipp Cannon's avatar
Kipp Cannon committed
712 713 714 715 716 717 718 719 720 721 722 723 724 725
		# point this to a LnLRDensity object containing the
		# zero-lag densities to mix zero-lag into the model.
		self.lnzerolagdensity = None

		# initialize a CoincRates object.  NOTE:  this is
		# potentially time-consuming.  the current implementation
		# includes hard-coded fast-paths for the standard
		# gstlal-based inspiral pipeline's coincidence and network
		# configurations, but if those change then doing this will
		# suck.  when scipy >= 0.19 becomes available on LDG
		# clusters this issue will go away (can use qhull's
		# algebraic geometry code for the probability
		# calculations).
		self.coinc_rates = snglcoinc.CoincRates(
726
			instruments = self.instruments,
727 728
			delta_t = self.delta_t,
			min_instruments = self.min_instruments
Kipp Cannon's avatar
Kipp Cannon committed
729 730 731 732 733 734
		)

	@property
	def segmentlists(self):
		return self.triggerrates.segmentlistdict()

735
	def __call__(self, segments, snrs, chi2s_over_snr2s, phase, dt, template_id):
Kipp Cannon's avatar
Kipp Cannon committed
736 737
		assert frozenset(segments) == self.instruments
		if len(snrs) < self.min_instruments:
738
			return NegInf
Kipp Cannon's avatar
Kipp Cannon committed
739

740 741 742 743 744 745
		# FIXME:  the +/-3600 s window thing is a temporary hack to
		# work around the problem of vetoes creating short segments
		# that have no triggers in them but that can have
		# injections recovered in them.  the +/- 3600 s window is
		# just a guess as to what might be sufficient to work
		# around it.  you might might to make this bigger.
Kipp Cannon's avatar
Kipp Cannon committed
746 747
		triggers_per_second_per_template = {}
		for instrument, seg in segments.items():
748
			triggers_per_second_per_template[instrument] = (self.triggerrates[instrument] & trigger_rate.ratebinlist([trigger_rate.ratebin(seg[1] - 3600., seg[1] + 3600., count = 0)])).density / len(self.template_ids)
Kipp Cannon's avatar
Kipp Cannon committed
749
		# sanity check rates
750
		assert all(triggers_per_second_per_template[instrument] for instrument in snrs), "impossible candidate in %s at %s when rates were %s triggers/s/template" % (", ".join(sorted(snrs)), ", ".join("%s s in %s" % (str(seg[1]), instrument) for instrument, seg in sorted(segments.items())), str(triggers_per_second_per_template))
Kipp Cannon's avatar
Kipp Cannon committed
751

Kipp Cannon's avatar
Kipp Cannon committed
752 753 754 755 756 757 758 759 760
		# P(t | noise) = (candidates per unit time @ t) / total
		# candidates.  by not normalizing by the total candidates
		# the return value can only ever be proportional to the
		# probability density, but we avoid the problem of the
		# ranking statistic definition changing on-the-fly while
		# running online, allowing candidates collected later to
		# have their ranking statistics compared meaningfully to
		# the values assigned to candidates collected earlier, when
		# the total number of candidates was smaller.
761
		lnP = math.log(sum(self.coinc_rates.strict_coinc_rates(**triggers_per_second_per_template).values()) * len(self.template_ids))
Kipp Cannon's avatar
Kipp Cannon committed
762 763 764 765 766

		# P(instruments | t, noise)
		lnP += self.coinc_rates.lnP_instruments(**triggers_per_second_per_template)[frozenset(snrs)]

		# evaluate dt and dphi parameters
767 768
		# NOTE: uniform and normalized so that the log should be zero, but there is no point in doing that
		# lnP += 0
Kipp Cannon's avatar
Kipp Cannon committed
769 770 771

		# evaluate the rest
		interps = self.interps
772
		return lnP + sum(interps["%s_snr_chi" % instrument](snrs[instrument], chi2_over_snr2) for instrument, chi2_over_snr2 in chi2s_over_snr2s.items())
Kipp Cannon's avatar
Kipp Cannon committed
773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799

	def __iadd__(self, other):
		super(LnNoiseDensity, self).__iadd__(other)
		self.triggerrates += other.triggerrates
		return self

	def copy(self):
		new = super(LnNoiseDensity, self).copy()
		new.triggerrates = self.triggerrates.copy()
		# NOTE:  lnzerolagdensity in the copy is reset to None by
		# this operation.  it is left as an exercise to the calling
		# code to re-connect it to the appropriate object if
		# desired.
		return new

	def mkinterps(self):
		#
		# override to mix zero-lag densities in if requested
		#

		if self.lnzerolagdensity is None:
			super(LnNoiseDensity, self).mkinterps()
		else:
			# same as parent class, but with .lnzerolagdensity
			# added
			self.interps = dict((key, (pdf + self.lnzerolagdensity.densities[key]).mkinterp()) for key, pdf in self.densities.items())

800
	def add_noise_model(self, number_of_events = 10000, prefactors_range = (2.0, 100.), df = 40, inv_snr_pow = 2.):
Kipp Cannon's avatar
Kipp Cannon committed
801 802 803 804 805 806 807
		#
		# populate snr,chi2 binnings with a slope to force
		# higher-SNR events to be assesed to be more significant
		# when in the regime beyond the edge of measured or even
		# extrapolated background.
		#

808 809 810 811 812 813 814 815 816 817 818
		# pick a canonical PDF to definine the binning (we assume
		# they're all the same and only compute this array once to
		# save time
		lnpdf = self.densities.values()[0]
		arr = numpy.zeros_like(lnpdf.array)

		snrindices, rcossindices = lnpdf.bins[self.snr_min:1e10, 1e-6:1e2]
		snr, dsnr = lnpdf.bins[0].centres()[snrindices], lnpdf.bins[0].upper()[snrindices] - lnpdf.bins[0].lower()[snrindices]
		rcoss, drcoss = lnpdf.bins[1].centres()[rcossindices], lnpdf.bins[1].upper()[rcossindices] - lnpdf.bins[1].lower()[rcossindices]

		prcoss = numpy.ones(len(rcoss))
819 820
		# This adds a faint power law that falls off just faster than GWs
		psnr = 1e-12 * snr**-6 #(1. + 10**6) / (1. + snr**6)
821
		psnr = numpy.outer(psnr, numpy.ones(len(rcoss)))
822 823
		# NOTE the magic numbers are just tuned from real data
		psnrdcoss = numpy.outer(numpy.exp(-4. * (snr - 4.5)**2) * dsnr, numpy.exp(-(rcoss - .06)**2 / (1e-4)) * drcoss)
824
		arr[snrindices, rcossindices] = psnrdcoss + psnr
825 826 827 828 829

		# normalize to the requested count.  give 99% of the
		# requested events to this portion of the model
		arr *= 0.99 * number_of_events / arr.sum()

Kipp Cannon's avatar
Kipp Cannon committed
830
		for lnpdf in self.densities.values():
831 832 833
			# add in the 99% noise model
			lnpdf.array += arr
			# add 1% from the "glitch model"
Kipp Cannon's avatar
Kipp Cannon committed
834 835 836 837 838 839 840 841 842 843
			inspiral_extrinsics.NumeratorSNRCHIPDF.add_signal_model(lnpdf, n = 0.01 * number_of_events, prefactors_range = prefactors_range, df = df, inv_snr_pow = inv_snr_pow, snr_min = self.snr_min)
			# re-normalize
			lnpdf.normalize()

	def candidate_count_model(self):
		"""
		Compute and return a prediction for the total number of
		noise candidates expected for each instrument
		combination.
		"""
844 845 846 847 848 849
		# assumes the trigger rate is uniformly distributed among
		# the templates and uniform in live time, calculates
		# coincidence rate assuming template exact-match
		# coincidence is required, then multiplies by the template
		# count to get total coincidence rate.
		return dict((instruments, count * len(self.template_ids)) for instruments, count in self.coinc_rates.marginalized_strict_coinc_counts(
Kipp Cannon's avatar
Kipp Cannon committed
850
			self.triggerrates.segmentlistdict(),
851
			**dict((instrument, rate / len(self.template_ids)) for instrument, rate in self.triggerrates.densities.items())
Kipp Cannon's avatar
Kipp Cannon committed
852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886
		).items())

	def random_params(self):
		"""
		Generator that yields an endless sequence of randomly
		generated candidate parameters.  NOTE: the parameters will
		be within the domain of the repsective binnings but are not
		drawn from the PDF stored in those binnings --- this is not
		an MCMC style sampler.  Each value in the sequence is a
		three-element tuple.  The first two elements of each tuple
		provide the *args and **kwargs values for calls to this PDF
		or the numerator PDF or the ranking statistic object.  The
		final is the natural logarithm (up to an arbitrary
		constant) of the PDF from which the parameters have been
		drawn evaluated at the point described by the *args and
		**kwargs.

		See also:

		random_sim_params()

		The sequence is suitable for input to the .ln_lr_samples()
		log likelihood ratio generator.
		"""
		snr_slope = 0.8 / len(self.instruments)**3

		snrchi2gens = dict((instrument, iter(self.densities["%s_snr_chi" % instrument].bins.randcoord(ns = (snr_slope, 1.), domain = (slice(self.snr_min, None), slice(None, None)))).next) for instrument in self.instruments)
		t_and_rate_gen = iter(self.triggerrates.random_uniform()).next
		t_offsets_gen = dict((instruments, self.coinc_rates.plausible_toas(instruments).next) for instruments in self.coinc_rates.all_instrument_combos)
		random_randint = random.randint
		random_sample = random.sample
		random_uniform = random.uniform
		segment = segments.segment
		twopi = 2. * math.pi
		ln_1_2 = math.log(0.5)
887 888
		lnP_template_id = -math.log(len(self.template_ids))
		template_ids = tuple(self.template_ids)
Kipp Cannon's avatar
Kipp Cannon committed
889 890 891 892 893 894 895 896 897 898 899 900 901 902
		def nCk(n, k):
			return math.factorial(n) // math.factorial(k) // math.factorial(n - k)
		while 1:
			t, rates, lnP_t = t_and_rate_gen()

			instruments = tuple(instrument for instrument, rate in rates.items() if rate > 0)
			if len(instruments) < self.min_instruments:
				# FIXME:  doing this invalidates lnP_t.  I
				# think the error is merely an overall
				# normalization error, though, and nothing
				# cares about the normalization
				continue
			# to pick instruments, we first pick an integer k
			# between m = min_instruments and n =
903
			# len(instruments) inclusively, then choose that
Kipp Cannon's avatar
Kipp Cannon committed
904 905 906 907 908 909
			# many unique names from among the available
			# instruments.  the probability of the outcome is
			#
			# = P(k) * P(selection | k)
			# = 1 / (n - m + 1) * 1 / nCk
			#
910 911
			# where nCk = number of k choices without
			# replacement from a collection of n things.
Kipp Cannon's avatar
Kipp Cannon committed
912 913 914 915 916
			k = random_randint(self.min_instruments, len(instruments))
			lnP_instruments = -math.log((len(instruments) - self.min_instruments + 1) * nCk(len(instruments), k))
			instruments = frozenset(random_sample(instruments, k))

			seq = sum((snrchi2gens[instrument]() for instrument in instruments), ())
917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934
			kwargs = {
				# FIXME: waveform duration hard-coded to
				# 10 s, generalize
				"segments": dict.fromkeys(self.instruments, segment(t - 10.0, t)),
				"snrs": dict((instrument, value[0]) for instrument, value in zip(instruments, seq[0::2])),
				"chi2s_over_snr2s": dict((instrument, value[1]) for instrument, value in zip(instruments, seq[0::2])),
				"phase": dict((instrument, random_uniform(0., twopi)) for instrument in instruments),
				# FIXME:  add dt to segments?  not
				# self-consistent if we don't, but doing so
				# screws up the test that was done to check
				# which instruments are on and off at "t"
				"dt": t_offsets_gen[instruments](),
				# FIXME random_params needs to be given a
				# meaningful template_id, but for now it is
				# not used in the likelihood-ratio
				# assignment so we don't care
				"template_id": random.choice(template_ids)
			}
Kipp Cannon's avatar
Kipp Cannon committed
935 936 937 938 939 940 941
			# NOTE:  I think the result of this sum is, in
			# fact, correctly normalized, but nothing requires
			# it to be (only that it be correct up to an
			# unknown constant) and I've not checked that it is
			# so the documentation doesn't promise that it is.
			# FIXME:  no, it's not normalized until the dt_dphi
			# bit is corrected for other than H1L1
942
			yield (), kwargs, sum(seq[1::2], lnP_t + lnP_instruments + lnP_template_id)
Kipp Cannon's avatar
Kipp Cannon committed
943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966

	def to_xml(self, name):
		xml = super(LnNoiseDensity, self).to_xml(name)
		xml.appendChild(self.triggerrates.to_xml(u"triggerrates"))
		return xml

	@classmethod
	def from_xml(cls, xml, name):
		xml = cls.get_xml_root(xml, name)
		self = super(LnNoiseDensity, cls).from_xml(xml, name)
		self.triggerrates = trigger_rate.triggerrates.from_xml(xml, u"triggerrates")
		self.triggerrates.coalesce()	# just in case
		return self


class DatalessLnNoiseDensity(LnNoiseDensity):
	"""
	Stripped-down version of LnNoiseDensity for use in estimating
	ranking statistics when no data has been collected from an
	instrument with which to define the ranking statistic.

	Used, for example, to implement low-significance candidate culls,
	etc.

967 968 969 970
	Assumes all available instruments are on and have the same horizon
	distance, and assess candidates based only on SNR and \chi^2
	distributions.
	"""
971 972 973 974 975 976 977 978

	DEFAULT_FILENAME = os.path.join(gstlal_config_paths["pkgdatadir"], "inspiral_datalesslndensity.xml.gz")

	@ligolw_array.use_in
	@ligolw_param.use_in
	class LIGOLWContentHandler(ligolw.LIGOLWContentHandler):
		pass

Kipp Cannon's avatar
Kipp Cannon committed
979 980
	def __init__(self, *args, **kwargs):
		super(DatalessLnNoiseDensity, self).__init__(*args, **kwargs)
981 982 983 984 985 986 987

		# install SNR, chi^2 PDF (one for all instruments)
		# FIXME:  make mass dependent
		self.densities = {
			"snr_chi": rate.BinnedLnPDF.from_xml(ligolw_utils.load_filename(self.DEFAULT_FILENAME, contenthandler = self.LIGOLWContentHandler), u"datalesslnnoisedensity")
		}

988

989
	def __call__(self, segments, snrs, chi2s_over_snr2s, phase, dt, template_id):
990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007
		# assume all instruments are on, 1 trigger per second per
		# template
		triggers_per_second_per_template = dict.fromkeys(segments, 1.)

		# P(t | noise) = (candidates per unit time @ t) / total
		# candidates.  by not normalizing by the total candidates
		# the return value can only ever be proportional to the
		# probability density, but we avoid the problem of the
		# ranking statistic definition changing on-the-fly while
		# running online, allowing candidates collected later to
		# have their ranking statistics compared meaningfully to
		# the values assigned to candidates collected earlier, when
		# the total number of candidates was smaller.
		lnP = math.log(sum(self.coinc_rates.strict_coinc_rates(**triggers_per_second_per_template).values()) * len(self.template_ids))

		# P(instruments | t, noise)
		lnP += self.coinc_rates.lnP_instruments(**triggers_per_second_per_template)[frozenset(snrs)]

1008 1009
		# evalute the (snr, \chi^2 | snr) PDFs (same for all
		# instruments)
1010
		interp = self.interps["snr_chi"]
1011
		return lnP + sum(interp(snrs[instrument], chi2_over_snr2) for instrument, chi2_over_snr2 in chi2s_over_snr2s.items())
Kipp Cannon's avatar
Kipp Cannon committed
1012 1013 1014 1015 1016

	def random_params(self):
		# won't work
		raise NotImplementedError

1017 1018 1019 1020 1021 1022 1023 1024 1025
	def __iadd__(self, other):
		raise NotImplementedError

	def increment(self, *args, **kwargs):
		raise NotImplementedError

	def copy(self):
		raise NotImplementedError

Kipp Cannon's avatar
Kipp Cannon committed
1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037
	def to_xml(self, name):
		# I/O not permitted:  the on-disk version would be
		# indistinguishable from a real ranking statistic and could
		# lead to accidents
		raise NotImplementedError

	@classmethod
	def from_xml(cls, xml, name):
		# I/O not permitted:  the on-disk version would be
		# indistinguishable from a real ranking statistic and could
		# lead to accidents
		raise NotImplementedError
1038 1039


Kipp Cannon's avatar
Kipp Cannon committed
1040
class OnlineFrankensteinLnNoiseDensity(LnNoiseDensity):
1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054
	"""
	Version of LnNoiseDensity with trigger rate data spliced in from
	another instance.  Used to solve a chicken-or-egg problem and
	assign ranking statistic values in an aonline anlysis.  NOTE:  the
	trigger rate data is not copied from the donor, instances of this
	class hold a reference to the donor's data, so as it is modified
	those modifications are immediately reflected here.

	For safety's sake, instances cannot be written to or read from
	files, cannot be marginalized together with other instances, nor
	accept updates from new data.
	"""
	@classmethod
	def splice(cls, src, rates_donor):
Kipp Cannon's avatar
Kipp Cannon committed
1055
		self = cls(src.template_ids, src.instruments, src.delta_t, min_instruments = src.min_instruments)
1056
		for key, lnpdf in src.densities.items():
1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083
			self.densities[key] = lnpdf.copy()
		# NOTE:  not a copy.  we hold a reference to the donor's
		# data so that as it is updated, we get the updates.
		self.triggerrates = rates_donor.triggerrates
		return self

	def __iadd__(self, other):
		raise NotImplementedError

	def increment(self, *args, **kwargs):
		raise NotImplementedError

	def copy(self):
		raise NotImplementedError

	def to_xml(self, name):
		# I/O not permitted:  the on-disk version would be
		# indistinguishable from a real ranking statistic and could
		# lead to accidents
		raise NotImplementedError

	@classmethod
	def from_xml(cls, xml, name):
		# I/O not permitted:  the on-disk version would be
		# indistinguishable from a real ranking statistic and could
		# lead to accidents
		raise NotImplementedError