diff git a/lalburst/configure.ac b/lalburst/configure.ac
index 3aff588b3cd0819a6a798e71be90e987678c2fa1..a87e298dceb4df868e82b24549d077204ae14edb 100644
 a/lalburst/configure.ac
+++ b/lalburst/configure.ac
@@ 72,10 +72,12 @@ MIN_LAL_VERSION="7.2.2"
MIN_LALMETAIO_VERSION="3.0.0"
MIN_LALSIMULATION_VERSION="4.0.0"
MIN_PYTHON_LIGO_LW_VERSION="1.7.0"
+MIN_SCIPY_VERSION="0.19"
AC_SUBST([MIN_LAL_VERSION])
AC_SUBST([MIN_LALMETAIO_VERSION])
AC_SUBST([MIN_LALSIMULATION_VERSION])
AC_SUBST([MIN_PYTHON_LIGO_LW_VERSION])
+AC_SUBST([MIN_SCIPY_VERSION])
AC_CANONICAL_HOST
diff git a/lalburst/debian/control.in b/lalburst/debian/control.in
index 3831a3d3bfa6b00dffa68cab04fa77f18f54b0d7..c8776c66e62539b67f81496896647fa08fb8d7fb 100644
 a/lalburst/debian/control.in
+++ b/lalburst/debian/control.in
@@ 24,7 +24,7 @@ BuildDepends:
python3numpy,
python3pil,
python3pytest,
 python3scipy,
+ python3scipy (>= @MIN_SCIPY_VERSION@~),
python3tqdm,
rsync,
swig (>= @MIN_SWIG_VERSION@)  swig3.0 (>= @MIN_SWIG_VERSION@),
@@ 88,7 +88,7 @@ Depends:
python3ligolw,
python3ligosegments,
python3matplotlib,
 python3scipy,
+ python3scipy (>= @MIN_SCIPY_VERSION@~),
python3tqdm,
Description: Python 3 bindings foar LALBurst
The LSC Algorithm Burst Library for gravitational wave data analysis.
@@ 108,7 +108,7 @@ Depends:
python3matplotlib,
python3numpy,
python3pil,
 python3scipy,
+ python3scipy (>= @MIN_SCIPY_VERSION@~),
python3tqdm,
Description: LSC Algorithm Library Burst
The LSC Algorithm Burst Library for gravitational wave data analysis.
diff git a/lalburst/lalburst.spec.in b/lalburst/lalburst.spec.in
index 270cb3ff3e1589eb8cb96792d629ddae234988f1..6c38f36ee8732c6fb9e9e995b0a93938b620b2ca 100644
 a/lalburst/lalburst.spec.in
+++ b/lalburst/lalburst.spec.in
@@ 47,7 +47,7 @@ BuildRequires: python%{python3_pkgversion}matplotlib
BuildRequires: python%{python3_pkgversion}numpy >= @MIN_NUMPY_VERSION@
BuildRequires: python%{python3_pkgversion}pillow
BuildRequires: python%{python3_pkgversion}pytest
BuildRequires: python%{python3_pkgversion}scipy
+BuildRequires: python%{python3_pkgversion}scipy >= @MIN_SCIPY_VERSION@
BuildRequires: python3rpmmacros
# octave
diff git a/lalburst/python/lalburst/snglcoinc.py b/lalburst/python/lalburst/snglcoinc.py
index 6930864df665930c8e090cf51b1d8c7a92ce8651..7fead29528ec057a1ea368cefc4c1c3ca6e95a76 100644
 a/lalburst/python/lalburst/snglcoinc.py
+++ b/lalburst/python/lalburst/snglcoinc.py
@@ 1606,82 +1606,112 @@ class CoincRates(object):
self.rate_factors = {}
for instruments in self.all_instrument_combos:
 # choose the instrument whose TOA forms the "epoch" of the
 # coinc. to improve the convergence rate this should be
 # the instrument with the smallest Cartesian product of
 # coincidence windows with other instruments (so that
 # coincidence with this instrument provides the tightest
 # prior constraint on the time differences between the
 # other instruments).
+ # choose the instrument whose event time forms the "epoch"
+ # of the coinc. which we choose is irrelevant
+ # frozenset of names
key = instruments
 anchor = min(instruments, key = lambda a: sum(math.log(self.tau[frozenset((a, b))]) for b in instruments  set([a])))
 instruments = tuple(instruments  set([anchor]))
 # the computation of a coincidence rate starts by computing
 # \mu_{1} * \mu_{2} ... \mu_{N} * 2 * \tau_{12} * 2 *
 # \tau_{13} ... 2 * \tau_{1N}. this is the rate at which
 # events from instrument 1 are coincident with events from
 # all of instruments 2...N. the factors of 2 are because
 # to be coincident the time difference can be anywhere in
 # [tau, +tau], so the size of the coincidence window is 2
 # tau. removing the factor of
+ # name, list of names
+ anchor, *instruments = instruments
+ # to understand the calculation consider the following
+ # diagram. this depicts the interinstrument time delays
+ # for candidates in a threedetector network.
#
 # \prod_{i} \mu_{i}
+ # 3
+ # t3t2 = const ^
+ # /  t2t1 = const
+ #  /  
+ # +/++ t3t1 = const
+ # /XXXXXXXXXXXXXXX
+ # /XXXXXXXXXXXXXXXX
+ # /XXXXXXXXXXXXXXXX
+ # 1> 2
+ # XXXXXXXXXXXXXXXX/
+ # XXXXXXXXXXXXXXXX/
+ # XXXXXXXXXXXXXXX/
+ # ++/+ t3t1 = const
+ #   / 
+ # t2t1 = const  /
+ # t3t2 = const
#
 # leaves
+ # in this example the time of the event seen in detector 1
+ # (arbitrary) is represented by the origin. the x axis
+ # indicates the time difference between the event seen in
+ # detector 2 and the event seen in detector 1, while the y
+ # axis depicts the time differenc between the event seen in
+ # detector 3 and the event seen in detector 1. in this
+ # parameter space, the coincidence windows take the form of
+ # halfspaces: for example, the requirement that t2t1
+ # <= const becomes two t2t1=const halfspace boundaries,
+ # in this case vertical lines parallel to the yaxis.
+ # similarly, the requirement that t3t1 <= const becomes
+ # a pair of t3t1=const halfspace boundaries parallel to
+ # the xaxis. the requirement that t3t2 <= const can
+ # also be shown in this space, that requirement becomes a
+ # third pair of halfspace boundaries, this time at 45
+ # degrees to the axes.
#
 # \prod_{i} 2 \tau_{1i}.
+ # in this 2dimensional parameter space, all valid 3way
+ # mutual coincidences must fall within the shaded region
+ # marked with "X"s, meaning their (t2t1) and (t3t1)
+ # values must fall in that area. this is the intersection
+ # of the 6 halfspaces.
#
 # in the N1 dimensional space defined by the time
 # differences between each instrument and the anchor
 # instrument, the coincidence windows between instruments
 # define pairs of halfspace boundaries. for the
 # coincidence windows between each instrument and the
 # anchor instrument these boundaries are perpendicular to
 # coordinate axes, while for other pairs of instruments
 # the coincidence windows correspond to planes angled at 45
 # degrees in various orientations. altogether they define
 # a convex polyhedron containing the origin.
+ # all of this generalizes to more than 3 detectors,
+ # becoming an (N1)dimensional space of the \Delta t's,
+ # with the allowed volume being the intersection of the
+ # halfspaces defined by the coincidence windows. it's
+ # aparent that the final quantity we seek here is the
+ # volume of the convex polyhedron defined by all of the
+ # constraints. this can be computed using the qhull
+ # library's halfspace intersection implementation.
#
 # the product of taus, above, is the volume of the
 # rectangular polyhedron defined by the anchor instrument
 # constraints alone. it's aparent that the final quantity
 # we seek here is the volume of the convex polyhedron
 # defined by all of the constraints. this can be computed
 # using the qhull library's halfspace intersection
 # implementation.
 #
 # the halfspace instersection code assumes constraints of
 # the form
+ # the halfspace instersection code requires constraints be
+ # provided to it in the form
#
# A x + b <= 0,
#
# where A has size (n constraints x m dimensions), where
 # for N instruments m = N1. each coincidence window
+ # for N instruments m = N1. in this form, x is the vector
+ # of time differences between the events seen in each of
+ # the detectors and the event seen in the "anchor"
+ # instrument at the origin. each coincidence window
# between an instrument, i, and the anchor imposes two
# constraints of the form
#
# +/t_{i}  \tau_{1i} <= 0
#
 # for a total of 2*(N1) constraints. each coincidence
 # window between a pair of (nonanchor) instruments imposes
 # two constraints of the form
+ # giving 2*(N1) constraints. additionally, each
+ # coincidence window between a pair of nonanchor
+ # instruments imposes two constraints of the form
#
# +/(t_{i}  t_{j})  \tau_{ij} <= 0
#
 # for a total (N1)*(N2) constraints. altogether there
 # are
+ # for an additional (N1)*(N2) constraints. altogether
+ # there are
#
# n = (N1)^2 + (N1) = m * (m + 1)
#
# constraints
+ #
+ # the "rate factor" we compute below is the volume in units
+ # of time^(N1) of the convex polyhedron described above.
+ # this is the proportionality constant giving us the rate
+ # of Nway mutual coincidences
+ #
+ # coinc rate = (polyhedron volume) \\prod_{i} \\mu_{i}.
+ #
+ # the volume has units of time^(N1), the \\mu have units
+ # of 1/time, there are N of them, so as hoped the rate of
+ # coincidences has units of 1/time.
+
if not instruments:
# oneinstrument case, noop
self.rate_factors[key] = 1.
elif len(instruments) == 1:
# twoinstrument (1D) case, don't use qhull
 self.rate_factors[key] = 2. * self.tau[frozenset((anchor, instruments[0]))]
+ self.rate_factors[key] = 2. * self.tau[key]
else:
 # three and more instrument (2D and
+ # three and more instrument (2D and
# higher) case
dimensions = len(instruments) # anchor not included
halfspaces = numpy.zeros((dimensions * (dimensions + 1), dimensions + 1), dtype = "double")
@@ 1703,67 +1733,7 @@ class CoincRates(object):
# the origin is in the interior
interior = numpy.zeros((len(instruments),), dtype = "double")
# compute volume
 try:
 self.rate_factors[key] = spatial.ConvexHull(spatial.HalfspaceIntersection(halfspaces, interior).intersections).volume
 except AttributeError:
 # fallthrough to old version
 pass
 else:
 # it worked, continue
 continue

 # old stonethrowing version in case qhull
 # is not available. FIXME: remove when we
 # are sure it's not needed. for each
 # instrument 2...N, the interval within
 # which an event is coincident with
 # instrument 1
 windows = tuple((self.tau[frozenset((anchor, instrument))], +self.tau[frozenset((anchor, instrument))]) for instrument in instruments)
 # preassemble a sequence of instrument
 # index pairs and the maximum allowed
 # \Delta t between them to avoid doing the
 # work associated with assembling the
 # sequence inside a loop
 ijseq = tuple((i, j, self.tau[frozenset((instruments[i], instruments[j]))]) for (i, j) in itertools.combinations(range(len(instruments)), 2))
 # compute the numerator and denominator of
 # the fraction of events coincident with
 # the anchor instrument that are also
 # mutually coincident. this is done by
 # picking a vector of allowed \Delta ts and
 # testing them against the coincidence
 # windows. the loop's exit criterion is
 # arrived at as follows. after d trials,
 # the number of successful outcomes is a
 # binomiallydistributed RV with variance =
 # d p (1  p) <= d/4 where p is the
 # probability of a successful outcome. we
 # quit when the ratio of the bound on the
 # standard deviation of the number of
 # successful outcomes (\sqrt{d/4}) to the
 # actual number of successful outcomes (n)
 # falls below rel accuracy: \sqrt{d/4} / n
 # < rel accuracy, or
 #
 # \sqrt{d} < 2 * rel accuracy * n
 #
 # note that if the true probability is 0,
 # so that n=0 identically, then the loop
 # will never terminate; from the nature of
 # the problem we know 0= two_epsilon * n:
 dt = tuple(random_uniform(*window) for window in windows)
 if all(abs(dt[i]  dt[j]) <= maxdt for i, j, maxdt in ijseq):
 n += 1
 d += 1
 self.rate_factors[key] = float(n) / float(d)
 for instrument in instruments:
 self.rate_factors[key] *= 2. * self.tau[frozenset((anchor, instrument))]

+ self.rate_factors[key] = spatial.ConvexHull(spatial.HalfspaceIntersection(halfspaces, interior).intersections).volume
# done computing rate_factors
@@ 1956,6 +1926,14 @@ class CoincRates(object):
>>> x = iter(coincrates.random_instruments(H1 = 0.001, L1 = 0.002, V1 = 0.003))
>>> x.next() # doctest: +SKIP
(frozenset(['H1', 'L1']), 3.738788683913535)
+
+ NOTE: the generator yields instrument combinations drawn
+ uniformly from the set of allowed instrument combinations,
+ but this is not considered part of the definition of the
+ behaviour of this generator and this implementation detail
+ should not be relied upon by calling code. It is stated
+ here simply for clarity in case it helps improve
+ optimization choices elsewhere.
"""
# guaranteed nonempty
lnP_instruments = self.lnP_instruments(**rates)