NstarTools.py

Go to the documentation of this file.

# Copyright (C) 2017 Gregory Ashton
#
# 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.
"""
 
Tools to estimate Nstar and related functions
 
"""
from __future__ import division, absolute_import, print_function
 
import pprint
import itertools
import logging
import numpy as np
 
try:
    import lal
except ImportError:
    raise ImportError("SWIG wrappings of LAL cannot be imported")
 
try:
    import lalpulsar
except ImportError:
    raise ImportError("SWIG wrappings of LALPulsar cannot be imported")
 
 
def _extract_data_from_prior(prior):
    """Calculate the input data from the prior
 
    Parameters
    ----------
    prior: dict
        For each key in 'Alpha', 'Delta', 'F0', 'F1', 'F2', either a scalar
        value (not searched over), or a pair of the upper and lower limit.
 
    Returns
    -------
    p : ndarray
        Matrix with columns being the edges of the uniform bounding box
    spindowns : int
        The number of spindowns
    sky : bool
        If true, search includes the sky position
    fiducial_freq : float
        Fidicual frequency
 
    """
    keys = ["Alpha", "Delta", "F0", "F1", "F2"]
    spindown_keys = keys[3:]
    sky_keys = keys[:2]
    lims = []
    lims_keys = []
    scalar_keys = []
    unfound_keys = []
    lims_idxs = []
    for i, key in enumerate(keys):
        array_like = isinstance(prior[key], (list, tuple, np.ndarray))
        if key not in keys:
            raise ValueError("Key {}={} not understood".format(key, prior[key]))
        elif array_like and len(prior[key]) == 2:
            lims.append([prior[key][0], prior[key][1]])
            lims_keys.append(key)
            lims_idxs.append(i)
        elif np.isscalar(prior[key]) or len(prior[key] == 1):
            scalar_keys.append(key)
            if key in sky_keys:
                lims.append([prior[key], 0])
        else:
            unfound_keys.append(key)
    logging.info(
        "Computing Nstar over dimensions {} while {} are scalar and {} unused".format(
            lims_keys, scalar_keys, unfound_keys
        )
    )
    lims = np.array(lims)
    lims_keys = np.array(lims_keys)
    base = lims[:, 0]
    p = [base]
    for i in lims_idxs:
        basex = base.copy()
        basex[i] = lims[i, 1]
        p.append(basex)
    spindowns = np.sum([np.sum(lims_keys == k) for k in spindown_keys])
    sky = any([key in lims_keys for key in sky_keys])
    fiducial_freq = np.max(prior["F0"])
 
    return np.array(p).T, spindowns, sky, fiducial_freq
 
 
def _convert_array_to_gsl_matrix(array):
    gsl_matrix = lal.gsl_matrix(*array.shape)
    gsl_matrix.data = array
    return gsl_matrix
 
 
def get_Nstar_estimate(
    nsegs,
    tref,
    minStartTime,
    maxStartTime,
    prior,
    detector_names,
    earth_ephem_file="earth00-40-DE405.dat.gz",
    sun_ephem_file="sun00-40-DE405.dat.gz",
):
    """Returns N* estimated from the super-sky metric
 
    Nstar is the approximate number of unit-mismatch templates, see
    https://dcc.ligo.org/P1700455 for further details.
 
    Parameters
    ----------
    nsegs : int
        Number of semi-coherent segments
    tref : int
        Reference time in GPS seconds
    minStartTime, maxStartTime : int
        Minimum and maximum SFT timestamps
    prior: dict
        For each key in 'Alpha', 'Delta', 'F0', 'F1', 'F2', either a scalar
        value (not searched over), or a pair of the upper and lower limit.
    detector_names : list
        List of detectors to average over, e.g. ['H1']
 
    Returns
    -------
    Nstar: int
        The estimated approximate number of templates to cover the prior
        parameter space at a mismatch of unity, assuming the normalised
        thickness is unity.
 
    Example
    --------
 
    An example for a directed search where the Nstar can be estimated from the
    metric directly. This estimate is used to define the frequency and
    spin-down uncertainty. The calculated estimate using `get_Nstar_estimate`
    agrees with this input value.
 
    >>> from lalpulsar import NstarTools
    >>> nsegs = 1
    >>> minStartTime = 1000000000
    >>> duration = 10 * 86400
    >>> maxStartTime = minStartTime + duration
    >>> tref = minStartTime + .5*duration
    >>> detector_names = ['H1']
    >>> F0 = 30
    >>> F1 = -1e-10
    >>> F2 = 0
    >>> Alpha = 0.5
    >>> Delta = 1.5
    >>> Nstar = 1e3
    >>> DeltaF0 = np.sqrt(Nstar) * np.sqrt(3)/(np.pi*duration)
    >>> DeltaF1 = np.sqrt(Nstar) * np.sqrt(180)/(np.pi*duration**2)
    >>> prior = {'F0': [F0-DeltaF0/2., F0+DeltaF0/2],
    >>>          'F1': [F1-DeltaF1/2., F1+DeltaF1/2],
    >>>          'F2': F2,
    >>>          'Alpha': Alpha
    >>>          'Delta': Delta}
    >>> print NstarTools.get_Nstar_estimate(
    >>>     nsegs, tref, minStartTime, maxStartTime, prior, detector_names)
    1000.00000009
 
    Note
    ----
 
    To see detailed information about each set of dimensions Nstar, add the
    following before calling get_Nstar_estimate()
    >>> import logging
    >>> logging.basicConfig(level=logging.DEBUG)
 
    """
    detector_names = list(detector_names)
 
    in_phys, spindowns, sky, fiducial_freq = _extract_data_from_prior(prior)
    out_rssky = np.zeros(in_phys.shape)
 
    # Convert to Alpha, Delta, Fn, Fn-1, Fn-2,.. ordering
    in_phys[2:] = in_phys[2:][::-1]
 
    in_phys = _convert_array_to_gsl_matrix(in_phys)
    out_rssky = _convert_array_to_gsl_matrix(out_rssky)
 
    tboundaries = np.linspace(minStartTime, maxStartTime, nsegs + 1)
 
    ref_time = lal.LIGOTimeGPS(tref)
    segments = lal.SegListCreate()
    for j in range(len(tboundaries) - 1):
        seg = lal.SegCreate(
            lal.LIGOTimeGPS(tboundaries[j]), lal.LIGOTimeGPS(tboundaries[j + 1]), j
        )
        lal.SegListAppend(segments, seg)
    detNames = lal.CreateStringVector(*detector_names)
    detectors = lalpulsar.MultiLALDetector()
    lalpulsar.ParseMultiLALDetector(detectors, detNames)
    detector_weights = None
    detector_motion = lalpulsar.DETMOTION_SPIN + lalpulsar.DETMOTION_ORBIT
    ephemeris = lalpulsar.InitBarycenter(earth_ephem_file, sun_ephem_file)
    try:
        SSkyMetric = lalpulsar.ComputeSuperskyMetrics(
            lalpulsar.SUPERSKY_METRIC_TYPE,
            spindowns,
            ref_time,
            segments,
            fiducial_freq,
            detectors,
            detector_weights,
            detector_motion,
            ephemeris,
        )
    except RuntimeError as e:
        logging.warning("Encountered run-time error {}".format(e))
        raise RuntimeError("Calculation of the SSkyMetric failed")
 
    if sky:
        i = 0
    else:
        i = 2
 
    lalpulsar.ConvertPhysicalToSuperskyPoints(
        out_rssky, in_phys, SSkyMetric.semi_rssky_transf
    )
 
    d = out_rssky.data
    g = SSkyMetric.semi_rssky_metric.data
 
    g[2:, 2:] = g[2:, 2:][::-1, ::-1]  # Convert to Alpha, Delta, F0, F1, F2
 
    # Remove sky if required
    g = g[i:, i:]
    d = d[i:, i:]
 
    # Calculate parallelepiped of differences
    parallelepiped = (d[:, 1:].T - d[:, 0]).T
 
    labels = np.array(["Alpha", "Delta", "F0", "F1", "F2"])[i:]
 
    # Compute Nstar over all allowed combinations then take maximum
    Nstars = {}
    for j in range(1, len(g) + 1):
        for idx in itertools.combinations(np.arange(len(g)), j):
            idx = list(idx)
            dV = np.abs(np.linalg.det(parallelepiped[idx][:, idx]))
            sqrtdetG = np.sqrt(np.abs(np.linalg.det(g[idx][:, idx])))
            Nstars["".join(labels[idx])] = sqrtdetG * dV
    logging.debug(pprint.pformat(Nstars))
    return np.max(Nstars.values())