Add new extremely fast direct monte-carlo method

bin/bank/pycbc_brute_bank

@@ -23,6 +23,7 @@ import numpy
 import h5py
 import logging
 import argparse
+import pickle
 import numpy.random
 from scipy.stats import gaussian_kde
 
@@ -46,12 +47,12 @@ parser.add_argument('--min',
  help='list of the minimum parameter values', nargs='+', type=float)
 parser.add_argument('--max',
  help='list of the maximum parameter values', nargs='+', type=float)
-parser.add_argument('--approximant', required=True,
- help='The waveform approximant to place')
+parser.add_argument('--approximant', required=False,
+ help='The waveform approximant to place.')
 parser.add_argument('--minimal-match', default=0.97, type=float)
 parser.add_argument('--buffer-length', default=2, type=float,
  help='size of waveform buffer in seconds')
-parser.add_argument('--max-signal-length', type= float, 
+parser.add_argument('--max-signal-length', type= float,
  help="When specified, it cuts the maximum length of the waveform model to the lengh provided")
 parser.add_argument('--sample-rate', default=2048, type=float,
  help='sample rate in seconds')
@@ -60,7 +61,7 @@ parser.add_argument('--enable-sigma-bound', action='store_true')
 parser.add_argument('--tau0-threshold', type=float)
 parser.add_argument('--permissive', action='store_true',
  help='Allow waveform generator to fail.')
-parser.add_argument('--placement-iterations', default=1000, type=int, 
+parser.add_argument('--placement-iterations', default=1000, type=int,
  help='Specify the number of attempts the bank should make when placing points. Use this option if the bank fails to place any points.')
 parser.add_argument('--seed', type=int, default=0)
 parser.add_argument('--tolerance', type=float)
@@ -69,6 +70,7 @@ parser.add_argument('--min-mchirp', type=float)
 parser.add_argument('--max-mchirp', type=float)
 parser.add_argument('--fixed-params', type=str, nargs='*')
 parser.add_argument('--fixed-values', type=float, nargs='*')
+parser.add_argument('--use-cross', action='store_true')
 parser.add_argument('--max-q', type=float)
 parser.add_argument('--tau0-crawl', type=float)
 parser.add_argument('--tau0-start', type=float)
@@ -279,7 +281,7 @@ class GenUniformWaveform(object):
  self.md = q._data[-100:]
  self.md2 = q._data[0:100]
 
- def generate(self, **kwds): 
+ def generate(self, **kwds):
  kwds.update(fdict)
  if args.max_signal_length is not None:
  flow = numpy.arange(self.f_lower, 100, .1)[::-1]
@@ -290,17 +292,19 @@ class GenUniformWaveform(object):
  f = flow[x]
  else:
  f = self.f_lower
+
+ if hasattr(kwds['approximant'], 'decode'):
+ kwds['approximant'] = kwds['approximant'].decode()
 
- kwds['f_lower'] = f
-
- if kwds['approximant'] in pycbc.waveform.fd_approximants(): 
+ if kwds['approximant'] in pycbc.waveform.fd_approximants():
  hp, hc = pycbc.waveform.get_fd_waveform(delta_f=self.delta_f,
  f_ref=10.0, **kwds)
-
-
+ if args.use_cross:
+ hp = hc
+
  if 'fratio' in kwds:
  hp = hc * kwds['fratio'] + hp * (1 - kwds['fratio'])
- 
+
  else:
  dt = 1.0 / args.sample_rate
  hp = pycbc.waveform.get_waveform_filter(
@@ -382,7 +386,8 @@ def draw(rtype):
  if args.input_config is not None and waveform_transforms is not None:
  params = transforms.apply_transforms(params, waveform_transforms)
 
- params['approximant'] = numpy.array([args.approximant]*size)
+ if args.approximant is not None:
+ params['approximant'] = numpy.array([args.approximant]*size)
 
  # Filter out stuff (kde method may also generate samples outside boundaries).
  l = None
@@ -445,7 +450,7 @@ def cdraw(rtype, ts, te):
  return None
 
  return p
- 
+
 tau0s = args.tau0_start
 tau0e = tau0s + args.tau0_crawl
 
@@ -471,7 +476,7 @@ while tau0s < args.tau0_end:
  if r > 10:
  conv = uconv
  kloop = 0
- 
+
  while ((kloop == 0) or (kconv / okconv) > .5) and len(bank) > 10:
  r += 1
  kloop += 1
@@ -480,8 +485,8 @@ while tau0s < args.tau0_end:
  bank, kconv = bank.check_params(gen, params, args.minimal_match)
  logging.info("%s: Round (K) (%s): %s Size: %s conv: %s added: %s",
  region, kloop, r, len(bank), kconv, len(bank) - blen)
- 
- 
+
+
  if uconv:
  logging.info('Ratio of convergences: %2.3f' % (kconv / (uconv)))
  logging.info('Progress: {:.0%} completed'.format(tau0e/args.tau0_end))
@@ -500,7 +505,7 @@ while tau0s < args.tau0_end:
  tau0e += args.tau0_crawl / 2
 
 o = h5py.File(args.output_file, 'w')
- 
+o.attrs['minimal_match'] = args.minimal_match
 for k in bank.keys():
  val = bank.key(k)
  if val.dtype.char == 'U':

pycbc/distributions/joint.py

@@ -266,12 +266,14 @@ def contains(self, params):
  params = self.apply_boundary_conditions(**params)
  result = True
  for dist in self.distributions:
- param_name = dist.params[0]
- contain_array = numpy.ones(len(params[param_name]), dtype=bool)
+ param_names = dist.params
+ vlen = len(params[param_names[0]])
+ contain_array = numpy.ones(vlen, dtype=bool)
  # note: enable `__contains__` in `pycbc.distributions.bounded`
  # to handle array-like input, it doesn't work now.
- for index, k in enumerate(params[param_name]):
- contain_array[index] = {param_name: k} in dist
+ for i in range(vlen):
+ data = {pname: params[pname][i] for pname in param_names}
+ contain_array[i] = data in dist
  result &= numpy.array(contain_array)
  result &= self.within_constraints(params)
  return result

pycbc/inference/io/base_hdf.py

@@ -613,8 +613,6 @@ def read_injections(self, group=None):
  group = '/'.join([group, self.injections_group])
  injset = InjectionSet(self.filename, hdf_group=group)
  injections = injset.table.view(FieldArray)
- # close the new open filehandler to self
- injset._injhandler.filehandler.close()
  return injections
 
  def write_command_line(self):

pycbc/inference/models/marginalized_gaussian_noise.py

@@ -19,6 +19,7 @@
 """
 
 import itertools
+import logging
 import numpy
 from scipy import special
 
@@ -207,6 +208,7 @@ class MarginalizedTime(DistMarg, BaseGaussianNoise):
  def __init__(self, variable_params,
  data, low_frequency_cutoff, psds=None,
  high_frequency_cutoff=None, normalize=False,
+ sample_rate=None,
  **kwargs):
 
  self.kwargs = kwargs
@@ -240,6 +242,13 @@ def __init__(self, variable_params,
  gates=self.gates, **kwargs['static_params'])
 
  self.dets = {}
+
+ if sample_rate is not None:
+ logging.info("Using %s sample rate for marginalization", sample_rate)
+ for det in self._whitened_data:
+ tlen = int(round(float(sample_rate) * self.whitened_data[det].duration))
+ flen = tlen // 2 + 1
+ self._whitened_data[det].resize(flen) 
 
  def _nowaveform_loglr(self):
  """Convenience function to set loglr values if no waveform generated.
@@ -325,15 +334,20 @@ def _loglr(self):
  for det in wfs:
  if det not in self.dets:
  self.dets[det] = Detector(det)
- fp, fc = self.dets[det].antenna_pattern(
- params['ra'],
- params['dec'],
- params['polarization'],
- params['tc'])
- dt = self.dets[det].time_delay_from_earth_center(params['ra'],
- params['dec'],
- params['tc'])
+
+ if self.precalc_antenna_factors:
+ fp, fc, dt = self.get_precalc_antenna_factors(det)
+ else:
+ fp, fc = self.dets[det].antenna_pattern(
+ params['ra'],
+ params['dec'],
+ params['polarization'],
+ params['tc'])
+ dt = self.dets[det].time_delay_from_earth_center(params['ra'],
+ params['dec'],
+ params['tc'])
  dtc = params['tc'] + dt
+
  cplx_hd = fp * cplx_hpd[det].at_time(dtc,
  interpolate='quadratic')
  cplx_hd += fc * cplx_hcd[det].at_time(dtc,

pycbc/inference/models/relbin.py

@@ -837,7 +837,8 @@ def _loglr(self):
  pol_phase = numpy.exp(-2.0j * p['polarization'])
 
  snrs = self.get_snr(wfs)
- self.snr_draw(snrs=snrs)
+ if not self.snr_draw(snrs=snrs):
+ return -numpy.inf
 
  for ifo in self.sh:
  if self.precalc_antenna_factors: