inspiral_lr.py 43.1 KB
 Kipp Cannon committed Nov 26, 2017 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 # # ============================================================================= #  Kipp Cannon committed Nov 06, 2018 29 30 31 32 33 try: from fpconst import NegInf except ImportError: # not all machines have fpconst installed NegInf = float("-inf")  Kipp Cannon committed Nov 26, 2017 34 35 import math import numpy  Kipp Cannon committed Jun 06, 2018 36 import os  Kipp Cannon committed Nov 26, 2017 37 38 39 40 import random from scipy import interpolate from scipy import stats import sys  Kipp Cannon committed Dec 27, 2017 41 import warnings  Chad Hanna committed Jan 19, 2019 42 import json  Kipp Cannon committed Nov 26, 2017 43 44   Kipp Cannon committed Feb 14, 2019 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  Duncan Meacher committed Jul 11, 2018 50 from ligo import segments  Kipp Cannon committed Nov 26, 2017 51 52 from gstlal.stats import horizonhistory from gstlal.stats import inspiral_extrinsics  Duncan Meacher committed Jul 16, 2018 53 from gstlal.stats import inspiral_intrinsics  Kipp Cannon committed Nov 26, 2017 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  Kipp Cannon committed Jun 06, 2018 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 committed Nov 26, 2017 68 69 __all__ = [ "LnSignalDensity",  Kipp Cannon committed Nov 06, 2018 70 71 72 73 74  "LnNoiseDensity", "DatalessLnSignalDensity", "DatalessLnNoiseDensity", "OnlineFrankensteinLnSignalDensity", "OnlineFrankensteinLnNoiseDensity"  Kipp Cannon committed Nov 26, 2017 75 76 77 ]  chad.hanna committed Oct 09, 2018 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 committed Nov 26, 2017 83 84 85 86 87 88 89 90 91 92 93 # # ============================================================================= # # Base Class # # ============================================================================= # class LnLRDensity(snglcoinc.LnLRDensity): # range of SNRs covered by this object  Chad Hanna committed Jul 07, 2019 94  snr_min = 4.0  Kipp Cannon committed Nov 26, 2017 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)))  Chad Hanna committed Oct 18, 2018 99  def __init__(self, template_ids, instruments, delta_t, min_instruments = 2):  Kipp Cannon committed Nov 26, 2017 100 101 102 103  # # check input #  Kipp Cannon committed Dec 06, 2017 104 105  if template_ids is not None and not template_ids: raise ValueError("template_ids cannot be empty")  Kipp Cannon committed Nov 26, 2017 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 #  Kipp Cannon committed Dec 06, 2017 117  self.template_ids = frozenset(template_ids) if template_ids is not None else template_ids  Kipp Cannon committed Nov 26, 2017 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)  Kipp Cannon committed Dec 06, 2017 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 committed Nov 26, 2017 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")  Kipp Cannon committed Dec 06, 2017 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 committed Nov 26, 2017 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):  Kipp Cannon committed Dec 06, 2017 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.  Kipp Cannon committed Sep 04, 2018 158  if event.template_id not in self.template_ids:  Kipp Cannon committed Dec 06, 2017 159  raise ValueError("event from wrong template")  Kipp Cannon committed Oct 01, 2018 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 committed Nov 26, 2017 181 182 183  self.densities["%s_snr_chi" % event.ifo].count[event.snr, event.chisq / event.snr**2.] += 1.0 def copy(self):  Chad Hanna committed Oct 18, 2018 184  new = type(self)(self.template_ids, self.instruments, self.delta_t, self.min_instruments)  Kipp Cannon committed Nov 26, 2017 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)  Kipp Cannon committed Oct 22, 2018 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))  Chad Hanna committed Oct 17, 2018 249  # FIXME this is not an ideal way to get only one into the file  Kipp Cannon committed Nov 26, 2017 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)  Kipp Cannon committed May 01, 2018 263  template_ids = ligolw_param.get_pyvalue(xml, u"template_ids")  Chad Hanna committed Oct 17, 2018 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.  Kipp Cannon committed May 01, 2018 271 272  if template_ids is not None: template_ids = frozenset(int(i) for i in template_ids.split(","))  Kipp Cannon committed Nov 26, 2017 273  self = cls(  Kipp Cannon committed May 01, 2018 274  template_ids = template_ids,  Kipp Cannon committed Nov 26, 2017 275 276  instruments = lsctables.instrumentsproperty.get(ligolw_param.get_pyvalue(xml, u"instruments")), delta_t = ligolw_param.get_pyvalue(xml, u"delta_t"),  Chad Hanna committed Oct 18, 2018 277  min_instruments = ligolw_param.get_pyvalue(xml, u"min_instruments")  Kipp Cannon committed Nov 26, 2017 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):  Kipp Cannon committed Dec 06, 2017 294  def __init__(self, *args, **kwargs):  Kipp Cannon committed Oct 12, 2018 295  population_model_file = kwargs.pop("population_model_file", None)  Ryan Michael Magee committed Jul 29, 2019 296  dtdphi_file = kwargs.pop("dtdphi_file", None)  Chad Hanna committed Jan 19, 2019 297  self.horizon_factors = kwargs.pop("horizon_factors", None)  Kipp Cannon committed Dec 06, 2017 298  super(LnSignalDensity, self).__init__(*args, **kwargs)  Kipp Cannon committed Nov 26, 2017 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  Kipp Cannon committed Dec 06, 2017 306  self.horizon_history = horizonhistory.HorizonHistories((instrument, horizonhistory.NearestLeafTree()) for instrument in self.instruments)  Kipp Cannon committed Oct 12, 2018 307  self.population_model_file = population_model_file  Ryan Michael Magee committed Jul 29, 2019 308  self.dtdphi_file = dtdphi_file  Chad Hanna committed Oct 17, 2018 309 310  # source population model # FIXME: introduce a mechanism for selecting the file  Chad Hanna committed Oct 18, 2018 311  self.population_model = inspiral_intrinsics.SourcePopulationModel(self.template_ids, filename = self.population_model_file)  Ryan Michael Magee committed Jul 29, 2019 312 313  if self.dtdphi_file is not None: self.InspiralExtrinsics = inspiral_extrinsics.InspiralExtrinsics(self.min_instruments, self.dtdphi_file)  Chad Hanna committed May 07, 2018 314   Chad Hanna committed Jan 19, 2019 315 316 317  def set_horizon_factors(self, horizon_factors): self.horizon_factors = horizon_factors  Kipp Cannon committed Sep 04, 2018 318  def __call__(self, segments, snrs, chi2s_over_snr2s, phase, dt, template_id):  Kipp Cannon committed Nov 26, 2017 319 320  assert frozenset(segments) == self.instruments if len(snrs) < self.min_instruments:  Kipp Cannon committed Nov 06, 2018 321  return NegInf  Kipp Cannon committed Nov 26, 2017 322   Kipp Cannon committed Dec 06, 2017 323 324  # use volume-weighted average horizon distance over # duration of event to estimate sensitivity  Kipp Cannon committed Nov 26, 2017 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())  chad.hanna committed Oct 09, 2018 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 committed Nov 26, 2017 334  # in lieu of properly normalizating this factor. we can't  chad.hanna committed Oct 09, 2018 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  Kipp Cannon committed Dec 06, 2017 348  if not horizon:  Kipp Cannon committed Nov 06, 2018 349  return NegInf  Chad Hanna committed Jan 19, 2019 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))  Kipp Cannon committed Dec 06, 2017 356   Chad Hanna committed May 07, 2018 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  Kipp Cannon committed Nov 06, 2018 362  return NegInf  Kipp Cannon committed Nov 26, 2017 363   Chad Hanna committed May 07, 2018 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:  Kipp Cannon committed Nov 06, 2018 369  return NegInf  Kipp Cannon committed Nov 26, 2017 370   Duncan Meacher committed Jul 16, 2018 371 372 373  # evaluate population model lnP += self.population_model.lnP_template_signal(template_id, max(snrs.values()))  Kipp Cannon committed Nov 27, 2017 374 375 376  # evalute the (snr, \chi^2 | snr) PDFs (same for all # instruments) interp = self.interps["snr_chi"]  Kipp Cannon committed Sep 04, 2018 377  return lnP + sum(interp(snrs[instrument], chi2_over_snr2) for instrument, chi2_over_snr2 in chi2s_over_snr2s.items())  Kipp Cannon committed Nov 26, 2017 378 379 380 381  def __iadd__(self, other): super(LnSignalDensity, self).__iadd__(other) self.horizon_history += other.horizon_history  Chad Hanna committed Oct 18, 2018 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  Ryan Michael Magee committed Jul 29, 2019 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  Chad Hanna committed Jan 19, 2019 390  if self.horizon_factors is not None and other.horizon_factors is not None and other.horizon_factors != self.horizon_factors:  Chad Hanna committed Jan 19, 2019 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")  Chad Hanna committed Jan 19, 2019 395 396  if self.horizon_factors is None and other.horizon_factors is not None: self.horizon_factors = other.horizon_factors  Chad Hanna committed Oct 18, 2018 397   Kipp Cannon committed Nov 26, 2017 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 committed Oct 12, 2018 406  new.population_model_file = self.population_model_file  Ryan Michael Magee committed Jul 29, 2019 407  new.dtdphi_file = self.dtdphi_file  Kipp Cannon committed Oct 12, 2018 408 409 410  # okay to use references because read-only data new.population_model = self.population_model new.InspiralExtrinsics = self.InspiralExtrinsics  Chad Hanna committed Jan 19, 2019 411  new.horizon_factors = self.horizon_factors  Kipp Cannon committed Nov 26, 2017 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())  Chad Hanna committed Aug 09, 2019 426  def add_signal_model(self, prefactors_range = (0.001, 0.30), df = 200, inv_snr_pow = 4.):  Kipp Cannon committed Nov 26, 2017 427 428  # normalize to 10 *mi*llion signals. this count makes the # density estimation code choose a suitable kernel size  Chad Hanna committed Aug 09, 2019 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 committed Nov 26, 2017 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.)  Kipp Cannon committed Dec 06, 2017 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 committed Nov 26, 2017 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  Kipp Cannon committed Feb 14, 2019 454  object (see ligo.lw.lsctables.SimInspiral for more  Kipp Cannon committed Nov 26, 2017 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  Kipp Cannon committed Sep 23, 2018 481  # FIXME need to add dt and dphi  Kipp Cannon committed Nov 26, 2017 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"))  Kipp Cannon committed Oct 22, 2018 554  xml.appendChild(ligolw_param.Param.from_pyvalue(u"population_model_file", self.population_model_file))  Chad Hanna committed Jan 19, 2019 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))  Ryan Michael Magee committed Jul 29, 2019 556  xml.appendChild(ligolw_param.Param.from_pyvalue(u"dtdphi_file", self.dtdphi_file))  Kipp Cannon committed Nov 26, 2017 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 committed Oct 26, 2018 564  self.population_model_file = ligolw_param.get_pyvalue(xml, u"population_model_file")  Ryan Michael Magee committed Jul 29, 2019 565  self.dtdphi_file = ligolw_param.get_pyvalue(xml, u"dtdphi_file")  Chad Hanna committed Jan 19, 2019 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 committed Oct 12, 2018 573  self.population_model = inspiral_intrinsics.SourcePopulationModel(self.template_ids, filename = self.population_model_file)  Ryan Michael Magee committed Jul 29, 2019 574 575  if self.dtdphi_file is not None: self.InspiralExtrinsics = inspiral_extrinsics.InspiralExtrinsics(self.min_instruments, self.dtdphi_file)  Kipp Cannon committed Nov 26, 2017 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.  Kipp Cannon committed Dec 06, 2017 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 committed Nov 26, 2017 591  """  Kipp Cannon committed Dec 06, 2017 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 committed Nov 26, 2017 596   Kipp Cannon committed Sep 04, 2018 597  def __call__(self, segments, snrs, chi2s_over_snr2s, phase, dt, template_id):  Kipp Cannon committed Dec 06, 2017 598 599  # evaluate P(t) \propto number of templates lnP = math.log(len(self.template_ids))  Kipp Cannon committed Nov 26, 2017 600   Kipp Cannon committed Nov 06, 2018 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 committed Nov 06, 2018 608 609  # raises ValueError when a needed horizon distance # is zero  Kipp Cannon committed Nov 06, 2018 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:  Kipp Cannon committed Nov 06, 2018 617  return NegInf  Kipp Cannon committed Nov 26, 2017 618   Duncan Meacher committed Jul 16, 2018 619 620 621  # evaluate population model lnP += self.population_model.lnP_template_signal(template_id, max(snrs.values()))  Kipp Cannon committed Dec 06, 2017 622 623 624  # evalute the (snr, \chi^2 | snr) PDFs (same for all # instruments) interp = self.interps["snr_chi"]  Kipp Cannon committed Sep 04, 2018 625  return lnP + sum(interp(snrs[instrument], chi2_over_snr2) for instrument, chi2_over_snr2 in chi2s_over_snr2s.items())  Kipp Cannon committed Nov 26, 2017 626 627 628 629  def __iadd__(self, other): raise NotImplementedError  Kipp Cannon committed Dec 06, 2017 630 631 632  def increment(self, *args, **kwargs): raise NotImplementedError  Kipp Cannon committed Nov 26, 2017 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 committed Nov 05, 2018 650 class OnlineFrankensteinLnSignalDensity(LnSignalDensity):  Kipp Cannon committed Apr 18, 2018 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):  chad.hanna committed Jan 19, 2019 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)  Cody Messick committed May 13, 2018 666  for key, lnpdf in src.densities.items():  Kipp Cannon committed Apr 18, 2018 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 committed Nov 26, 2017 696 697 698 699 700 701 702 703 704 705 # # ============================================================================= # # Denominator Density # # ============================================================================= # class LnNoiseDensity(LnLRDensity):  Kipp Cannon committed Dec 06, 2017 706 707  def __init__(self, *args, **kwargs): super(LnNoiseDensity, self).__init__(*args, **kwargs)  Kipp Cannon committed Nov 26, 2017 708 709 710  # record of trigger counts vs time for all instruments in # the network  Kipp Cannon committed Dec 06, 2017 711  self.triggerrates = trigger_rate.triggerrates((instrument, trigger_rate.ratebinlist()) for instrument in self.instruments)  Kipp Cannon committed Nov 26, 2017 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(  Chad Hanna committed May 07, 2018 726  instruments = self.instruments,  Kipp Cannon committed Dec 06, 2017 727 728  delta_t = self.delta_t, min_instruments = self.min_instruments  Kipp Cannon committed Nov 26, 2017 729 730 731 732 733 734  ) @property def segmentlists(self): return self.triggerrates.segmentlistdict()  Kipp Cannon committed Sep 04, 2018 735  def __call__(self, segments, snrs, chi2s_over_snr2s, phase, dt, template_id):  Kipp Cannon committed Nov 26, 2017 736 737  assert frozenset(segments) == self.instruments if len(snrs) < self.min_instruments:  Kipp Cannon committed Nov 06, 2018 738  return NegInf  Kipp Cannon committed Nov 26, 2017 739   Kipp Cannon committed Jan 31, 2018 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 committed Nov 26, 2017 746 747  triggers_per_second_per_template = {} for instrument, seg in segments.items():  CHAD RICHARD HANNA committed Jul 02, 2019 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 committed Nov 26, 2017 749  # sanity check rates  CHAD RICHARD HANNA committed Jul 02, 2019 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 committed Nov 26, 2017 751   Kipp Cannon committed Dec 04, 2017 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.  Kipp Cannon committed Dec 06, 2017 761  lnP = math.log(sum(self.coinc_rates.strict_coinc_rates(**triggers_per_second_per_template).values()) * len(self.template_ids))  Kipp Cannon committed Nov 26, 2017 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  chad.hanna committed Nov 17, 2018 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 committed Nov 26, 2017 769 770 771  # evaluate the rest interps = self.interps  Kipp Cannon committed Sep 04, 2018 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 committed Nov 26, 2017 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())  Chad Hanna committed Jun 14, 2019 800  def add_noise_model(self, number_of_events = 10000, prefactors_range = (2.0, 100.), df = 40, inv_snr_pow = 2.):  Kipp Cannon committed Nov 26, 2017 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. #  Kipp Cannon committed Dec 06, 2017 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))  Chad Hanna committed Jun 14, 2019 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)  Chad Hanna committed Jun 06, 2019 821  psnr = numpy.outer(psnr, numpy.ones(len(rcoss)))  Chad Hanna committed Jun 14, 2019 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)  Chad Hanna committed Jun 06, 2019 824  arr[snrindices, rcossindices] = psnrdcoss + psnr  Kipp Cannon committed Dec 06, 2017 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 committed Nov 26, 2017 830  for lnpdf in self.densities.values():  Kipp Cannon committed Dec 06, 2017 831 832 833  # add in the 99% noise model lnpdf.array += arr # add 1% from the "glitch model"  Kipp Cannon committed Nov 26, 2017 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. """  Kipp Cannon committed Dec 06, 2017 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 committed Nov 26, 2017 850  self.triggerrates.segmentlistdict(),  Kipp Cannon committed Dec 06, 2017 851  **dict((instrument, rate / len(self.template_ids)) for instrument, rate in self.triggerrates.densities.items())  Kipp Cannon committed Nov 26, 2017 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)  Kipp Cannon committed May 18, 2018 887 888  lnP_template_id = -math.log(len(self.template_ids)) template_ids = tuple(self.template_ids)  Kipp Cannon committed Nov 26, 2017 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 =  Kipp Cannon committed Sep 04, 2018 903  # len(instruments) inclusively, then choose that  Kipp Cannon committed Nov 26, 2017 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 #  Kipp Cannon committed Sep 04, 2018 910 911  # where nCk = number of k choices without # replacement from a collection of n things.  Kipp Cannon committed Nov 26, 2017 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), ())  Kipp Cannon committed Sep 04, 2018 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 committed Nov 26, 2017 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  Kipp Cannon committed May 18, 2018 942  yield (), kwargs, sum(seq[1::2], lnP_t + lnP_instruments + lnP_template_id)  Kipp Cannon committed Nov 26, 2017 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.  Kipp Cannon committed Dec 06, 2017 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. """  Kipp Cannon committed Jun 06, 2018 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 committed Nov 26, 2017 979 980  def __init__(self, *args, **kwargs): super(DatalessLnNoiseDensity, self).__init__(*args, **kwargs)  Kipp Cannon committed Jun 06, 2018 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") }  Kipp Cannon committed Dec 06, 2017 988   Kipp Cannon committed Sep 04, 2018 989  def __call__(self, segments, snrs, chi2s_over_snr2s, phase, dt, template_id):  Kipp Cannon committed Dec 06, 2017 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)]  Kipp Cannon committed Sep 04, 2018 1008 1009  # evalute the (snr, \chi^2 | snr) PDFs (same for all # instruments)  Kipp Cannon committed Jun 06, 2018 1010  interp = self.interps["snr_chi"]  Kipp Cannon committed Sep 04, 2018 1011  return lnP + sum(interp(snrs[instrument], chi2_over_snr2) for instrument, chi2_over_snr2 in chi2s_over_snr2s.items())  Kipp Cannon committed Nov 26, 2017 1012 1013 1014 1015 1016  def random_params(self): # won't work raise NotImplementedError  Kipp Cannon committed Dec 06, 2017 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 committed Nov 26, 2017 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  Kipp Cannon committed Apr 18, 2018 1038 1039   Kipp Cannon committed Nov 05, 2018 1040 class OnlineFrankensteinLnNoiseDensity(LnNoiseDensity):  Kipp Cannon committed Apr 18, 2018 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 committed Nov 05, 2018 1055  self = cls(src.template_ids, src.instruments, src.delta_t, min_instruments = src.min_instruments)  Cody Messick committed May 13, 2018 1056  for key, lnpdf in src.densities.items():  Kipp Cannon committed Apr 18, 2018 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