core.py 48.7 KB
Newer Older
1001
        self.tCWFstatMapVersion = 'lal'
1002
        self.cudaDeviceName = None
Gregory Ashton's avatar
Gregory Ashton committed
1003
1004
1005
1006
1007
1008
1009
        self.init_computefstatistic_single_point()
        self.init_semicoherent_parameters()

    def init_semicoherent_parameters(self):
        logging.info(('Initialising semicoherent parameters from {} to {} in'
                      ' {} segments').format(
            self.minStartTime, self.maxStartTime, self.nsegs))
David Keitel's avatar
David Keitel committed
1010
        self.transientWindowType = 'rect'
Gregory Ashton's avatar
Gregory Ashton committed
1011
1012
1013
        self.whatToCompute = lalpulsar.FSTATQ_2F+lalpulsar.FSTATQ_ATOMS_PER_DET
        self.tboundaries = np.linspace(self.minStartTime, self.maxStartTime,
                                       self.nsegs+1)
1014
        self.Tcoh = self.tboundaries[1] - self.tboundaries[0]
Gregory Ashton's avatar
Gregory Ashton committed
1015

1016
1017
1018
1019
1020
1021
1022
1023
1024
        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
1025

1026
    def get_semicoherent_twoF(
1027
1028
1029
1030
1031
            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
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
        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
1049
        #if not self.transientWindowType:
1050
1051
1052
1053
1054
1055
1056
1057
        #    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
1058
1059

        detStat = 0
1060
1061
        if record_segments:
            self.detStat_per_segment = []
Gregory Ashton's avatar
Gregory Ashton committed
1062

1063
1064
1065
        self.windowRange.tau = int(self.Tcoh)  # TYPE UINT4
        for tstart in self.tboundaries[:-1]:
            d_detStat = self._get_per_segment_det_stat(tstart)
1066
1067
1068
            detStat += d_detStat
            if record_segments:
                self.detStat_per_segment.append(d_detStat)
Gregory Ashton's avatar
Gregory Ashton committed
1069
1070
1071

        return detStat

1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
    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
1101

1102
class SemiCoherentGlitchSearch(ComputeFstat):
Gregory Ashton's avatar
Gregory Ashton committed
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
    """ 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,
1113
                 nglitch=1, sftfilepattern=None, theta0_idx=0, BSGL=False,
1114
                 minCoverFreq=None, maxCoverFreq=None, assumeSqrtSX=None,
1115
                 detectors=None, SSBprec=None, injectSources=None):
Gregory Ashton's avatar
Gregory Ashton committed
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
        """
        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).
1126
1127
1128
        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
1129
1130
1131
1132
1133
1134
1135
1136
1137
        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)
1138
        self.set_ephemeris_files()
David Keitel's avatar
David Keitel committed
1139
1140
1141
        self.transientWindowType = 'rect'
        self.t0Band  = None
        self.tauBand = None
1142
        self.tCWFstatMapVersion = 'lal'
1143
        self.cudaDeviceName = None
David Keitel's avatar
David Keitel committed
1144
        self.binary  = False
Gregory Ashton's avatar
Gregory Ashton committed
1145
1146
        self.init_computefstatistic_single_point()

1147
    def get_semicoherent_nglitch_twoF(self, F0, F1, F2, Alpha, Delta, *args):
Gregory Ashton's avatar
Gregory Ashton committed
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
        """ 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)

1161
        thetas = self._calculate_thetas(theta, delta_thetas, tboundaries,
1162
                                        theta0_idx=self.theta0_idx)
Gregory Ashton's avatar
Gregory Ashton committed
1163
1164
1165
1166

        twoFSum = 0
        for i, theta_i_at_tref in enumerate(thetas):
            ts, te = tboundaries[i], tboundaries[i+1]
1167
            if te - ts > 1800:
1168
1169
1170
1171
                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
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188

        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

1189
        theta_at_glitch = self._shift_coefficients(theta, tglitch - tref)
Gregory Ashton's avatar
Gregory Ashton committed
1190
        theta_post_glitch_at_glitch = theta_at_glitch + delta_theta
1191
        theta_post_glitch = self._shift_coefficients(
Gregory Ashton's avatar
Gregory Ashton committed
1192
1193
            theta_post_glitch_at_glitch, tref - tglitch)

1194
        twoFsegA = self.get_fullycoherent_twoF(
Gregory Ashton's avatar
Gregory Ashton committed
1195
1196
1197
1198
1199
1200
            self.minStartTime, tglitch, theta[0], theta[1], theta[2], Alpha,
            Delta)

        if tglitch == self.maxStartTime:
            return twoFsegA

1201
        twoFsegB = self.get_fullycoherent_twoF(
Gregory Ashton's avatar
Gregory Ashton committed
1202
1203
1204
1205
1206
            tglitch, self.maxStartTime, theta_post_glitch[0],
            theta_post_glitch[1], theta_post_glitch[2], Alpha,
            Delta)

        return twoFsegA + twoFsegB