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 @@ Build-Depends: python3-numpy, python3-pil, python3-pytest, - python3-scipy, + python3-scipy (>= @MIN_SCIPY_VERSION@~), python3-tqdm, rsync, swig (>= @MIN_SWIG_VERSION@) | swig3.0 (>= @MIN_SWIG_VERSION@), @@ -88,7 +88,7 @@ Depends: python3-ligo-lw, python3-ligo-segments, python3-matplotlib, - python3-scipy, + python3-scipy (>= @MIN_SCIPY_VERSION@~), python3-tqdm, Description: Python 3 bindings foar LALBurst The LSC Algorithm Burst Library for gravitational wave data analysis. @@ -108,7 +108,7 @@ Depends: python3-matplotlib, python3-numpy, python3-pil, - python3-scipy, + python3-scipy (>= @MIN_SCIPY_VERSION@~), python3-tqdm, 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: python3-rpm-macros # 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 inter-instrument time delays + # for candidates in a three-detector network. # - # \prod_{i} \mu_{i} + # 3 + # t3-t2 = const ^ + # / | t2-t1 = const + # | / | | + # --+-/------+--------+-- t3-t1 = const + # |/XXXXXXX|XXXXXXXX| + # /XXXXXXXX|XXXXXXXX| + # /|XXXXXXXX|XXXXXXXX| + # ---|--------1--------|---> 2 + # |XXXXXXXX|XXXXXXXX|/ + # |XXXXXXXX|XXXXXXXX/ + # |XXXXXXXX|XXXXXXX/| + # --+--------+------/-+-- t3-t1 = const + # | | / | + # t2-t1 = const | / + # t3-t2 = 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 + # half-spaces: for example, the requirement that |t2-t1| + # <= const becomes two t2-t1=const half-space boundaries, + # in this case vertical lines parallel to the y-axis. + # similarly, the requirement that |t3-t1| <= const becomes + # a pair of t3-t1=const half-space boundaries parallel to + # the x-axis. the requirement that |t3-t2| <= const can + # also be shown in this space, that requirement becomes a + # third pair of half-space boundaries, this time at 45 + # degrees to the axes. # - # \prod_{i} 2 \tau_{1i}. + # in this 2-dimensional parameter space, all valid 3-way + # mutual coincidences must fall within the shaded region + # marked with "X"s, meaning their (t2-t1) and (t3-t1) + # values must fall in that area. this is the intersection + # of the 6 half-spaces. # - # in the N-1 dimensional space defined by the time - # differences between each instrument and the anchor - # instrument, the coincidence windows between instruments - # define pairs of half-space boundaries. for the - # coincidence windows between each instrument and the - # anchor instrument these boundaries are perpendicular to - # co-ordinate 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 (N-1)-dimensional space of the \Delta t's, + # with the allowed volume being the intersection of the + # half-spaces 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 half-space 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 half-space intersection - # implementation. - # - # the half-space instersection code assumes constraints of - # the form + # the half-space 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 = N-1. each coincidence window + # for N instruments m = N-1. 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*(N-1) constraints. each coincidence - # window between a pair of (non-anchor) instruments imposes - # two constraints of the form + # giving 2*(N-1) constraints. additionally, each + # coincidence window between a pair of non-anchor + # instruments imposes two constraints of the form # # +/-(t_{i} - t_{j}) - \tau_{ij} <= 0 # - # for a total (N-1)*(N-2) constraints. altogether there - # are + # for an additional (N-1)*(N-2) constraints. altogether + # there are # # n = (N-1)^2 + (N-1) = m * (m + 1) # # constraints + # + # the "rate factor" we compute below is the volume in units + # of time^(N-1) of the convex polyhedron described above. + # this is the proportionality constant giving us the rate + # of N-way mutual coincidences + # + # coinc rate = (polyhedron volume) \\prod_{i} \\mu_{i}. + # + # the volume has units of time^(N-1), 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: # one-instrument case, no-op self.rate_factors[key] = 1. elif len(instruments) == 1: # two-instrument (1-D) 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 (2-D and + # three and more instrument (2-D 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: - # fall-through to old version - pass - else: - # it worked, continue - continue - - # old stone-throwing 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) - # pre-assemble 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 - # binomially-distributed 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 non-empty lnP_instruments = self.lnP_instruments(**rates)