Python PALutils.dailyAverage方法代码示例

# 需要导入模块: import PALutils [as 别名]
# 或者: from PALutils import dailyAverage [as 别名]
    def addpulsar(self, parfile, timfile, DMOFF=None, DMXOFF=None, dailyAverage=False):

        """
        Add another pulsar to the HDF5 file, given a tempo2 par and tim file.

        @param parfile: tempo2 par file
        @param timfile: tempo2 tim file
        @param DMOFF: Option to turn off DMMODEL fitting
        @param DMOFF: Option to turn off DMX fitting
        @param dailyAverage: Option to perform daily averaging to reduce the number
                             of points by consructing daily averaged TOAs that have
                             one residual per day per frequency band. (This has only
                             been tested on NANOGrav data thus far.)

        """

        # Check whether the two files exist
        if not os.path.isfile(parfile) or not os.path.isfile(timfile):
            raise IOError, "Cannot find parfile (%s) or timfile (%s)!" % (parfile, timfile)
        assert(self.filename != None), "ERROR: HDF5 file not set!"

        # 'a' means: read/write if exists, create otherwise
        self.h5file = h5.File(self.filename, 'a')

        if "Model" in self.h5file:
            self.h5file.close()
            self.h5file = None
            raise IOError, "model already available in '%s'. Refusing to add data" % (self.filename)

        # Create the data subgroup if it does not exist
        if "Data" in self.h5file:
            datagroup = self.h5file["Data"]
        else:
            datagroup = self.h5file.create_group("Data")

        # Load pulsar data from the JPL Cython tempo2 library
        t2pulsar = t2.tempopulsar(parfile, timfile)

        # do multiple fits
        #t2pulsar.fit(iters=10)

        # turn off DMMODEL fitting
        if DMOFF is not None:
            t2pulsar['DMMODEL'].fit = False
        
        # turn off DMX fitting
        if DMXOFF is not None:
            DMXFlag = False
            print 'Turning off DMX fitting and turning DM fitting on'
            for par in t2pulsar.pars:
                if 'DMX' in par:
                    print par
                    t2pulsar[par].fit = False
                    t2pulsar['DM'].fit = True
                    DMXFlag = True
            if DMXFlag== False: 
                print 'NO DMX for pulsar {0}'.format(t2pulsar.name)

        # refit 5 times to make sure we are converged
        t2pulsar.fit(iters=5)

        # Create the pulsar subgroup if it does not exist
        if "Pulsars" in datagroup:
            pulsarsgroup = datagroup["Pulsars"]
        else:
            pulsarsgroup = datagroup.create_group("Pulsars")

        # Look up the name of the pulsar, and see if it exist
        if t2pulsar.name in pulsarsgroup:
            self.h5file.close()
            raise IOError, "%s already exists in %s!" % (t2pulsar.name, self.filename)

        pulsarsgroup = pulsarsgroup.create_group(t2pulsar.name)

        # Read the data from the tempo2 structure.
        designmatrix = np.double(t2pulsar.designmatrix())
        residuals = np.double(t2pulsar.residuals())    
        toas = np.double(t2pulsar.toas())
        errs = np.double(t2pulsar.toaerrs*1e-6)
        pname = t2pulsar.name

        try:    # if tim file has frequencies
            freqs = np.double(t2pulsar.freqs)
        except AttributeError: 
            freqs = 0

        try:    # if tim file has frequency band flags
            bands = t2pulsar.flags['B']
        except KeyError:
            bands = 0

        # if doing daily averaging
        if dailyAverage:

            # get average quantities
            toas, qmatrix, errs, dmatrix, freqs, bands = PALutils.dailyAverage(t2pulsar)

            # construct new daily averaged residuals and designmatrix
            residuals = np.dot(qmatrix, residuals)
            designmatrix = np.dot(qmatrix, dmatrix)
#.........这里部分代码省略.........

# 需要导入模块: import PALutils [as 别名]
# 或者: from PALutils import dailyAverage [as 别名]
    def addInverseCovFromNoiseFile(self, parfile, timfile, noisefile, DMOFF=None, DMXOFF=None, dailyAverage=False):
        """
        
        Add noise covariance matrix after timing model subtraction.

        """

        # Check whether the two files exist
        if not os.path.isfile(parfile) or not os.path.isfile(timfile):
            raise IOError, "Cannot find parfile (%s) or timfile (%s)!" % (parfile, timfile)
        assert(self.filename != None), "ERROR: HDF5 file not set!"

        # 'a' means: read/write if exists, create otherwise
        self.h5file = h5.File(self.filename, 'a')

        # Create the data subgroup if it does not exist
        if "Data" in self.h5file:
            datagroup = self.h5file["Data"]
        else:
            raise IOError, "Cannot add noise parameters if Data group does not exist!"

        # Load pulsar data from the JPL Cython tempo2 library
        t2pulsar = t2.tempopulsar(parfile, timfile)
        
        # turn off DMMODEL fitting
        if DMOFF is not None:
            t2pulsar['DMMODEL'].fit = False

        # turn off DMX fitting
        if DMXOFF is not None:
            DMXFlag = False
            print 'Turning off DMX fitting and turning DM fitting on'
            for par in t2pulsar.pars:
                if 'DMX' in par:
                    t2pulsar[par].fit = False
                    t2pulsar['DM'].fit = True
                    DMXFlag = True
            if DMXFlag== False: 
                print 'NO DMX for pulsar {0}'.format(t2pulsar.name)

        # refit 5 times to make sure we are converged
        t2pulsar.fit(iters=5)

        # Create the pulsar subgroup if it does not exist
        if "Pulsars" in datagroup:
            pulsarsgroup = datagroup["Pulsars"]
        else:
            raise IOError, "Cannot add noise parameters if pulsar group does not exist!"

        # Look up the name of the pulsar, and see if it exist
        if t2pulsar.name in pulsarsgroup:
            pass
        else:
            raise IOError, "%s must already exists in %s to add noise parameters!"\
                    % (t2pulsar.name, self.filename)

        pulsarsgroup = pulsarsgroup[t2pulsar.name]

        # first create G matrix from design matrix and toas
        designmatrix = np.double(t2pulsar.designmatrix())
        toas = np.double(t2pulsar.toas()*86400)
        errs = np.double(t2pulsar.toaerrs*1e-6)

        # if doing daily averaging
        if dailyAverage:

            # get average quantities
            toas, qmatrix, errs, dmatrix, freqs, bands = PALutils.dailyAverage(t2pulsar)

            # construct new daily averaged residuals and designmatrix
            toas *= 86400
            designmatrix = np.dot(qmatrix, dmatrix)
        
        G = PALutils.createGmatrix(designmatrix)

        # create matrix of time lags
        tm = PALutils.createTimeLags(toas, toas, round=True)

        # now read noise file to get model and parameters
        file = open(noisefile,'r')

        fH = None
        tau = None
        DMAmp = None
        DMgam = None
 
        for line in file.readlines():
            
            # default parameters for different models other than pure PL
            key = line.split()[0]

            # get amplitude
            if "Amp" == key:
                Amp = float(line.split()[-1])

            # get spectral index
            elif "gam" == key:
                gam = float(line.split()[-1])
            
            # get efac
#.........这里部分代码省略.........