optimal_setup_functions.py

"""

Provides functions to aid in calculating the optimal setup based on the metric
volume estimates.

"""

import logging
import numpy as np
import scipy.optimize
import lal
import lalpulsar


def get_optimal_setup(
        R, Nsegs0, tref, minStartTime, maxStartTime, DeltaOmega,
        DeltaFs, fiducial_freq, detector_names, earth_ephem, sun_ephem):
    logging.info('Calculating optimal setup for R={}, Nsegs0={}'.format(
        R, Nsegs0))

    V_0 = get_V_estimate(
        Nsegs0, tref, minStartTime, maxStartTime, DeltaOmega, DeltaFs,
        fiducial_freq, detector_names, earth_ephem, sun_ephem)
    logging.info('Stage {}, nsegs={}, V={}'.format(0, Nsegs0, V_0))

    nsegs_vals = [Nsegs0]
    V_vals = [V_0]

    i = 0
    nsegs_i = Nsegs0
    while nsegs_i > 1:
        nsegs_i, V_i = get_nsegs_ip1(
            nsegs_i, R, tref, minStartTime, maxStartTime, DeltaOmega,
            DeltaFs, fiducial_freq, detector_names, earth_ephem, sun_ephem)
        nsegs_vals.append(nsegs_i)
        V_vals.append(V_i)
        i += 1
        logging.info(
            'Stage {}, nsegs={}, V={}'.format(i, nsegs_i, V_i))

    return nsegs_vals, V_vals


def get_nsegs_ip1(
        nsegs_i, R, tref, minStartTime, maxStartTime, DeltaOmega,
        DeltaFs, fiducial_freq, detector_names, earth_ephem, sun_ephem):

    log10R = np.log10(R)
    log10Vi = np.log10(get_V_estimate(
        nsegs_i, tref, minStartTime, maxStartTime, DeltaOmega, DeltaFs,
        fiducial_freq, detector_names, earth_ephem, sun_ephem))

    def f(nsegs_ip1):
        if nsegs_ip1[0] > nsegs_i:
            return 1e6
        if nsegs_ip1[0] < 0:
            return 1e6
        nsegs_ip1 = int(nsegs_ip1[0])
        if nsegs_ip1 == 0:
            nsegs_ip1 = 1
        Vip1 = get_V_estimate(
            nsegs_ip1, tref, minStartTime, maxStartTime, DeltaOmega,
            DeltaFs, fiducial_freq, detector_names, earth_ephem, sun_ephem)
        if Vip1[0] is None:
            return 1e6
        else:
            log10Vip1 = np.log10(Vip1)
            return np.abs(log10Vi[0] + log10R - log10Vip1[0])
    res = scipy.optimize.minimize(f, .5*nsegs_i, method='Powell', tol=0.1,
                                  options={'maxiter': 10})
    nsegs_ip1 = int(res.x)
    if nsegs_ip1 == 0:
        nsegs_ip1 = 1
    if res.success:
        return nsegs_ip1, get_V_estimate(
            nsegs_ip1, tref, minStartTime, maxStartTime, DeltaOmega, DeltaFs,
            fiducial_freq, detector_names, earth_ephem, sun_ephem)
    else:
        raise ValueError('Optimisation unsuccesful')


def get_V_estimate(
        nsegs, tref, minStartTime, maxStartTime, DeltaOmega, DeltaFs,
        fiducial_freq, detector_names, earth_ephem, sun_ephem):
    """ Returns V, Vsky, Vpe estimated from the super-sky metric

    Parameters
    ----------
    nsegs: int
        Number of semi-coherent segments
    tref: int
        Reference time in GPS seconds
    minStartTime, maxStartTime: int
        Minimum and maximum SFT timestamps
    DeltaOmega: float
        Solid angle of the sky-patch
    DeltaFs: array
        Array of [DeltaF0, DeltaF1, ...], length determines the number of
        spin-down terms.
    fiducial_freq: float
        Fidicual frequency
    detector_names: array
        Array of detectors to average over
    earth_ephem, sun_ephem: st
        Paths to the ephemeris files

    """
    spindowns = len(DeltaFs) - 1
    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, sun_ephem)
    try:
        SSkyMetric = lalpulsar.ComputeSuperskyMetrics(
            spindowns, ref_time, segments, fiducial_freq, detectors,
            detector_weights, detector_motion, ephemeris)
    except RuntimeError as e:
        logging.debug('Encountered run-time error {}'.format(e))
        return None, None, None

    sqrtdetG_SKY = np.sqrt(np.linalg.det(
        SSkyMetric.semi_rssky_metric.data[:2, :2]))
    sqrtdetG_PE = np.sqrt(np.linalg.det(
        SSkyMetric.semi_rssky_metric.data[2:, 2:]))

    Vsky = .5*sqrtdetG_SKY*DeltaOmega
    Vpe = sqrtdetG_PE * np.prod(DeltaFs)
    if Vsky == 0:
        Vsky = 1
    if Vpe == 0:
        Vpe = 1
    return (Vsky * Vpe, Vsky, Vpe)