core.py 48.1 KB
Newer Older
Gregory Ashton's avatar
Gregory Ashton committed
1001
1002
        self.tboundaries = np.linspace(self.minStartTime, self.maxStartTime,
                                       self.nsegs+1)
1003
        self.Tcoh = self.tboundaries[1] - self.tboundaries[0]
Gregory Ashton's avatar
Gregory Ashton committed
1004

1005
1006
1007
1008
1009
1010
1011
1012
1013
        if hasattr(self, 'SFT_timestamps'):
            if self.tboundaries[0] < self.SFT_timestamps[0]:
                logging.debug(
                    'Semi-coherent start time {} before first SFT timestamp {}'
                    .format(self.tboundaries[0], self.SFT_timestamps[0]))
            if self.tboundaries[-1] > self.SFT_timestamps[-1]:
                logging.debug(
                    'Semi-coherent end time {} after last SFT timestamp {}'
                    .format(self.tboundaries[-1], self.SFT_timestamps[-1]))
Gregory Ashton's avatar
Gregory Ashton committed
1014

1015
    def get_semicoherent_twoF(
1016
1017
1018
1019
1020
            self, F0, F1, F2, Alpha, Delta, asini=None,
            period=None, ecc=None, tp=None, argp=None,
            record_segments=False):
        """ Returns twoF or ln(BSGL) semi-coherently at a single point """

Gregory Ashton's avatar
Gregory Ashton committed
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
        self.PulsarDopplerParams.fkdot = np.array([F0, F1, F2, 0, 0, 0, 0])
        self.PulsarDopplerParams.Alpha = Alpha
        self.PulsarDopplerParams.Delta = Delta
        if self.binary:
            self.PulsarDopplerParams.asini = asini
            self.PulsarDopplerParams.period = period
            self.PulsarDopplerParams.ecc = ecc
            self.PulsarDopplerParams.tp = tp
            self.PulsarDopplerParams.argp = argp

        lalpulsar.ComputeFstat(self.FstatResults,
                               self.FstatInput,
                               self.PulsarDopplerParams,
                               1,
                               self.whatToCompute
                               )

David Keitel's avatar
David Keitel committed
1038
        #if not self.transientWindowType:
1039
1040
1041
1042
1043
1044
1045
1046
        #    if self.BSGL is False:
        #        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
1047
1048

        detStat = 0
1049
1050
        if record_segments:
            self.detStat_per_segment = []
Gregory Ashton's avatar
Gregory Ashton committed
1051

1052
1053
1054
        self.windowRange.tau = int(self.Tcoh)  # TYPE UINT4
        for tstart in self.tboundaries[:-1]:
            d_detStat = self._get_per_segment_det_stat(tstart)
1055
1056
1057
            detStat += d_detStat
            if record_segments:
                self.detStat_per_segment.append(d_detStat)
Gregory Ashton's avatar
Gregory Ashton committed
1058
1059
1060

        return detStat

1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
    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
1090

1091
class SemiCoherentGlitchSearch(ComputeFstat):
Gregory Ashton's avatar
Gregory Ashton committed
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
    """ 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,
1102
                 nglitch=1, sftfilepattern=None, theta0_idx=0, BSGL=False,
1103
                 minCoverFreq=None, maxCoverFreq=None, assumeSqrtSX=None,
1104
                 detectors=None, SSBprec=None, injectSources=None):
Gregory Ashton's avatar
Gregory Ashton committed
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
        """
        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).
1115
1116
1117
        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
1118
1119
1120
1121
1122
1123
1124
1125
1126
        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)
1127
        self.set_ephemeris_files()
David Keitel's avatar
David Keitel committed
1128
1129
1130
        self.transientWindowType = 'rect'
        self.t0Band  = None
        self.tauBand = None
1131
        self.tCWFstatMapVersion = 'lal'
1132
        self.cudaDeviceName = None
David Keitel's avatar
David Keitel committed
1133
        self.binary  = False
Gregory Ashton's avatar
Gregory Ashton committed
1134
1135
        self.init_computefstatistic_single_point()

1136
    def get_semicoherent_nglitch_twoF(self, F0, F1, F2, Alpha, Delta, *args):
Gregory Ashton's avatar
Gregory Ashton committed
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
        """ 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)

1150
        thetas = self._calculate_thetas(theta, delta_thetas, tboundaries,
1151
                                        theta0_idx=self.theta0_idx)
Gregory Ashton's avatar
Gregory Ashton committed
1152
1153
1154
1155

        twoFSum = 0
        for i, theta_i_at_tref in enumerate(thetas):
            ts, te = tboundaries[i], tboundaries[i+1]
1156
            if te - ts > 1800:
1157
1158
1159
1160
                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
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177

        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

1178
        theta_at_glitch = self._shift_coefficients(theta, tglitch - tref)
Gregory Ashton's avatar
Gregory Ashton committed
1179
        theta_post_glitch_at_glitch = theta_at_glitch + delta_theta
1180
        theta_post_glitch = self._shift_coefficients(
Gregory Ashton's avatar
Gregory Ashton committed
1181
1182
            theta_post_glitch_at_glitch, tref - tglitch)

1183
        twoFsegA = self.get_fullycoherent_twoF(
Gregory Ashton's avatar
Gregory Ashton committed
1184
1185
1186
1187
1188
1189
            self.minStartTime, tglitch, theta[0], theta[1], theta[2], Alpha,
            Delta)

        if tglitch == self.maxStartTime:
            return twoFsegA

1190
        twoFsegB = self.get_fullycoherent_twoF(
Gregory Ashton's avatar
Gregory Ashton committed
1191
1192
1193
1194
1195
            tglitch, self.maxStartTime, theta_post_glitch[0],
            theta_post_glitch[1], theta_post_glitch[2], Alpha,
            Delta)

        return twoFsegA + twoFsegB
For faster browsing, not all history is shown. View entire blame