core.py 46.2 KB
Newer Older
1001
1002
1003
1004
1005
1006
1007
        #        return self.FstatResults.twoF[0]
        #    twoF = np.float(self.FstatResults.twoF[0])
        #    self.twoFX[0] = self.FstatResults.twoFPerDet(0)
        #    self.twoFX[1] = self.FstatResults.twoFPerDet(1)
        #    log10_BSGL = lalpulsar.ComputeBSGL(twoF, self.twoFX,
        #                                       self.BSGLSetup)
        #    return log10_BSGL/np.log10(np.exp(1))
Gregory Ashton's avatar
Gregory Ashton committed
1008
1009

        detStat = 0
1010
1011
        if record_segments:
            self.detStat_per_segment = []
Gregory Ashton's avatar
Gregory Ashton committed
1012

1013
1014
1015
        self.windowRange.tau = int(self.Tcoh)  # TYPE UINT4
        for tstart in self.tboundaries[:-1]:
            d_detStat = self._get_per_segment_det_stat(tstart)
1016
1017
1018
            detStat += d_detStat
            if record_segments:
                self.detStat_per_segment.append(d_detStat)
Gregory Ashton's avatar
Gregory Ashton committed
1019
1020
1021

        return detStat

1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
    def _get_per_segment_det_stat(self, tstart):
        self.windowRange.t0 = int(tstart)  # TYPE UINT4

        FS = lalpulsar.ComputeTransientFstatMap(
            self.FstatResults.multiFatoms[0], self.windowRange, False)

        if self.BSGL is False:
            d_detStat = 2*FS.F_mn.data[0][0]
        else:
            FstatResults_single = copy.copy(self.FstatResults)
            FstatResults_single.lenth = 1
            FstatResults_single.data = self.FstatResults.multiFatoms[0].data[0]
            FS0 = lalpulsar.ComputeTransientFstatMap(
                FstatResults_single.multiFatoms[0], self.windowRange, False)
            FstatResults_single.data = self.FstatResults.multiFatoms[0].data[1]
            FS1 = lalpulsar.ComputeTransientFstatMap(
                FstatResults_single.multiFatoms[0], self.windowRange, False)

            self.twoFX[0] = 2*FS0.F_mn.data[0][0]
            self.twoFX[1] = 2*FS1.F_mn.data[0][0]
            log10_BSGL = lalpulsar.ComputeBSGL(
                    2*FS.F_mn.data[0][0], self.twoFX, self.BSGLSetup)
            d_detStat = log10_BSGL/np.log10(np.exp(1))
        if np.isnan(d_detStat):
            logging.debug('NaNs in semi-coherent twoF treated as zero')
            d_detStat = 0

        return d_detStat

Gregory Ashton's avatar
Gregory Ashton committed
1051

1052
class SemiCoherentGlitchSearch(ComputeFstat):
Gregory Ashton's avatar
Gregory Ashton committed
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
    """ A semi-coherent glitch search

    This implements a basic `semi-coherent glitch F-stat in which the data
    is divided into segments either side of the proposed glitches and the
    fully-coherent F-stat in each segment is summed to give the semi-coherent
    F-stat
    """

    @helper_functions.initializer
    def __init__(self, label, outdir, tref, minStartTime, maxStartTime,
1063
                 nglitch=1, sftfilepattern=None, theta0_idx=0, BSGL=False,
1064
                 minCoverFreq=None, maxCoverFreq=None, assumeSqrtSX=None,
1065
                 detectors=None, SSBprec=None, injectSources=None):
Gregory Ashton's avatar
Gregory Ashton committed
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
        """
        Parameters
        ----------
        label, outdir: str
            A label and directory to read/write data from/to.
        tref, minStartTime, maxStartTime: int
            GPS seconds of the reference time, and start and end of the data.
        nglitch: int
            The (fixed) number of glitches; this can zero, but occasionally
            this causes issue (in which case just use ComputeFstat).
1076
1077
1078
        sftfilepattern: str
            Pattern to match SFTs using wildcards (*?) and ranges [0-9];
            mutiple patterns can be given separated by colons.
Gregory Ashton's avatar
Gregory Ashton committed
1079
1080
1081
1082
1083
1084
1085
1086
1087
        theta0_idx, int
            Index (zero-based) of which segment the theta refers to - uyseful
            if providing a tight prior on theta to allow the signal to jump
            too theta (and not just from)

        For all other parameters, see pyfstat.ComputeFStat.
        """

        self.fs_file_name = "{}/{}_FS.dat".format(self.outdir, self.label)
1088
        self.set_ephemeris_files()
David Keitel's avatar
David Keitel committed
1089
1090
1091
1092
        self.transientWindowType = 'rect'
        self.t0Band  = None
        self.tauBand = None
        self.binary  = False
Gregory Ashton's avatar
Gregory Ashton committed
1093
1094
        self.init_computefstatistic_single_point()

1095
    def get_semicoherent_nglitch_twoF(self, F0, F1, F2, Alpha, Delta, *args):
Gregory Ashton's avatar
Gregory Ashton committed
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
        """ Returns the semi-coherent glitch summed twoF """

        args = list(args)
        tboundaries = ([self.minStartTime] + args[-self.nglitch:]
                       + [self.maxStartTime])
        delta_F0s = args[-3*self.nglitch:-2*self.nglitch]
        delta_F1s = args[-2*self.nglitch:-self.nglitch]
        delta_F2 = np.zeros(len(delta_F0s))
        delta_phi = np.zeros(len(delta_F0s))
        theta = [0, F0, F1, F2]
        delta_thetas = np.atleast_2d(
                np.array([delta_phi, delta_F0s, delta_F1s, delta_F2]).T)

1109
        thetas = self._calculate_thetas(theta, delta_thetas, tboundaries,
1110
                                        theta0_idx=self.theta0_idx)
Gregory Ashton's avatar
Gregory Ashton committed
1111
1112
1113
1114

        twoFSum = 0
        for i, theta_i_at_tref in enumerate(thetas):
            ts, te = tboundaries[i], tboundaries[i+1]
1115
            if te - ts > 1800:
1116
1117
1118
1119
                twoFVal = self.get_fullycoherent_twoF(
                    ts, te, theta_i_at_tref[1], theta_i_at_tref[2],
                    theta_i_at_tref[3], Alpha, Delta)
                twoFSum += twoFVal
Gregory Ashton's avatar
Gregory Ashton committed
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136

        if np.isfinite(twoFSum):
            return twoFSum
        else:
            return -np.inf

    def compute_glitch_fstat_single(self, F0, F1, F2, Alpha, Delta, delta_F0,
                                    delta_F1, tglitch):
        """ Returns the semi-coherent glitch summed twoF for nglitch=1

        Note: OBSOLETE, used only for testing
        """

        theta = [F0, F1, F2]
        delta_theta = [delta_F0, delta_F1, 0]
        tref = self.tref

1137
        theta_at_glitch = self._shift_coefficients(theta, tglitch - tref)
Gregory Ashton's avatar
Gregory Ashton committed
1138
        theta_post_glitch_at_glitch = theta_at_glitch + delta_theta
1139
        theta_post_glitch = self._shift_coefficients(
Gregory Ashton's avatar
Gregory Ashton committed
1140
1141
            theta_post_glitch_at_glitch, tref - tglitch)

1142
        twoFsegA = self.get_fullycoherent_twoF(
Gregory Ashton's avatar
Gregory Ashton committed
1143
1144
1145
1146
1147
1148
            self.minStartTime, tglitch, theta[0], theta[1], theta[2], Alpha,
            Delta)

        if tglitch == self.maxStartTime:
            return twoFsegA

1149
        twoFsegB = self.get_fullycoherent_twoF(
Gregory Ashton's avatar
Gregory Ashton committed
1150
1151
1152
1153
1154
            tglitch, self.maxStartTime, theta_post_glitch[0],
            theta_post_glitch[1], theta_post_glitch[2], Alpha,
            Delta)

        return twoFsegA + twoFsegB