Python PALutils类代码示例

    def myloglike(cube, ndim, nparams):

        efac = cube[0]
        equad = 10**cube[1]
        gam = cube[2]
        A = 10**cube[3]

        fs = np.zeros(args.ss)
        for ii in range(args.ss):
            fs[ii] = 10**cube[4+ii]
        rho2 = np.zeros(args.ss)
        for ii in range(args.ss):
            rho2[ii] = cube[ii+4+args.ss]
        
        # check to make sure frequencies are ordered
        ordered = np.all([fs[ii] < fs[ii+1] for ii in range(args.ss-1)])

        if ordered:

            F1 = list(PALutils.createfourierdesignmatrix(psr.toas, args.nmodes).T)
            tmp, f = PALutils.createfourierdesignmatrix(psr.toas, args.nmodes, freq=True)
            for ii in range(args.ss):
                F1.append(np.cos(2*np.pi*fs[ii]*psr.toas))
                F1.append(np.sin(2*np.pi*fs[ii]*psr.toas))

            F = np.array(F1).T
            F = np.dot(proj, F)

            # compute rho from A and gam# compute total time span of data
            Tspan = psr.toas.max() - psr.toas.min()

            # get power spectrum coefficients
            f1yr = 1/3.16e7
            rho = list(np.log10(A**2/12/np.pi**2 * f1yr**(gam-3) * f**(-gam)/Tspan))

            # compute total rho
            for ii in range(args.ss):
                rho.append(rho2[ii])

            loglike = PALLikelihoods.lentatiMarginalizedLike(psr, F, s, np.array(rho), efac**2, equad)
            
        else:

            loglike = -np.inf

        #print efac, rho, loglike

        return loglike

def loglike(x):

    start = time.time()

    # parameter values
    theta = x[0]
    phi = x[1]
    freq = 10 ** x[2]
    mc = 10 ** x[3]
    dist = 10 ** x[4]
    psi = x[5]
    inc = x[6]
    phase = x[7]

    pdist = x[8:]

    # generate list of waveforms for all puslars
    s = PALutils.createResidualsFast(
        psr, theta, phi, mc, dist, freq, phase, psi, inc, pdist=pdist, evolve=False, phase_approx=True
    )

    loglike = 0
    for ct, p in enumerate(psr):

        diff = p.res - s[ct]
        loglike += -0.5 * (logdetTerm[ct] + np.dot(diff, np.dot(p.invCov, diff)))

    if np.isnan(loglike):
        print "NaN log-likelihood. Not good..."
        return -np.inf
    else:
        return loglike

def loglike(x):

    efac = np.ones(npsr)
    equad = np.zeros(npsr)
    theta = x[0]
    phi = x[1]
    f = 10**x[2]
    mc = 10**x[3]
    dist = 10**x[4]
    psi = x[5]
    inc = x[6]
    phase = x[7]
    pdist = x[8:(8+npsr)]
        
    loglike = 0
    for ct, p in enumerate(psr):

        # make waveform with pulsar term
        s = PALutils.createResiduals(p, theta, phi, mc, dist, f, phase, psi, inc, pdist=pdist[ct], \
                 psrTerm=True)

        # project onto white noise basis
        s = np.dot(projList[ct], s)

        loglike += np.sum(p.res*s/(efac[ct]*Diag[ct] + equad[ct]**2))
        loglike -= 0.5 * np.sum(s**2/(efac[ct]*Diag[ct] + equad[ct]**2))
   
    if np.isnan(loglike):
        print 'NaN log-likelihood. Not good...'
        return -np.inf
    else:
        return loglike

    def myloglike(cube, ndim, nparams):

        efac = cube[0]
        equad = 10**cube[1]

        fs = np.zeros(args.ss)
        for ii in range(args.ss):
            fs[ii] = 10**cube[2+ii]
        rho = np.zeros(args.nmodes+args.ss)
        for ii in range(args.nmodes+args.ss):
            rho[ii] = cube[ii+2+args.ss]

        F1 = list(PALutils.createfourierdesignmatrix(psr.toas, args.nmodes).T)
        for ii in range(args.ss):
            F1.append(np.cos(2*np.pi*fs[ii]*psr.toas))
            F1.append(np.sin(2*np.pi*fs[ii]*psr.toas))

        F = np.array(F1).T

        F = np.dot(proj, F)
       
        loglike = PALLikelihoods.lentatiMarginalizedLike(psr, F, s, rho, efac, equad)

        #print efac, rho, loglike

        return loglike

def optStat(psr, ORF, gam=4.33333):
    """
    Computes the Optimal statistic as defined in Chamberlin, Creighton, Demorest et al (2013)

    @param psr: List of pulsar object instances
    @param ORF: Vector of pairwise overlap reduction values
    @param gam: Power Spectral index of GBW (default = 13/3, ie SMBMBs)

    @return: Opt: Optimal statistic value (A_gw^2)
    @return: sigma: 1-sigma uncertanty on Optimal statistic
    @return: snr: signal-to-noise ratio of cross correlations

    """

    #TODO: maybe compute ORF in code instead of reading it in. Would be less
    # of a risk but a bit slower...

    k = 0
    npsr = len(psr)
    top = 0
    bot = 0
    for ll in xrange(0, npsr):
        for kk in xrange(ll+1, npsr):

            # form matrix of toa residuals and compute SigmaIJ
            tm = PALutils.createTimeLags(psr[ll].toas, psr[kk].toas)

            # create cross covariance matrix without overall amplitude A^2
            SIJ = ORF[k]/2 * PALutils.createRedNoiseCovarianceMatrix(tm, 1, gam)
            
            # construct numerator and denominator of optimal statistic
            bot += np.trace(np.dot(psr[ll].invCov, np.dot(SIJ, np.dot(psr[kk].invCov, SIJ.T))))
            top += np.dot(psr[ll].res, np.dot(psr[ll].invCov, np.dot(SIJ, \
                        np.dot(psr[kk].invCov, psr[kk].res))))
            k+=1

    # compute optimal statistic
    Opt = top/bot
    
    # compute uncertainty
    sigma = 1/np.sqrt(bot)

    # compute SNR
    snr = top/np.sqrt(bot)

    # return optimal statistic and snr
    return Opt, sigma, snr

def upperLimitFunc(h, fstat_ref, freq, nreal):
    """
    Compute the value of the fstat for a range of parameters, with fixed
    amplitude over many realizations.

    @param h: value of the strain amplitude to keep constant
    @param fstat_ref: value of fstat for real data set
    @param freq: GW frequency
    @param nreal: number of realizations

    """
    Tmaxyr = np.array([(p.toas.max() - p.toas.min())/3.16e7 for p in psr]).max()
    count = 0
    for ii in range(nreal):

        # draw parameter values
        gwtheta = np.arccos(np.random.uniform(-1, 1))
        gwphi = np.random.uniform(0, 2*np.pi)
        gwphase = np.random.uniform(0, 2*np.pi)
        gwinc = np.arccos(np.random.uniform(0, 1))
        gwpsi = np.random.uniform(-np.pi/4, np.pi/4)

        # check to make sure source has not coalesced during observation time
        coal = True
        while coal:
            gwmc = 10**np.random.uniform(7, 10)
            tcoal = 2e6 * (gwmc/1e8)**(-5/3) * (freq/1e-8)**(-8/3)
            if tcoal > Tmaxyr:
                coal = False

        # determine distance in order to keep strain fixed
        gwdist = 4 * np.sqrt(2/5) * (gwmc*4.9e-6)**(5/3) * (np.pi*freq)**(2/3) / h

        # convert back to Mpc
        gwdist /= 1.0267e14

        # create residuals 
        for ct,p in enumerate(psr):
            inducedRes = PALutils.createResiduals(p, gwtheta, gwphi, gwmc, gwdist, \
                            freq, gwphase, gwpsi, gwinc, evolve=True)
 
            # replace residuals in pulsar object
            p.res = res[ct] + np.dot(R[ct], inducedRes)

        # compute f-statistic
        fpstat = PALLikelihoods.fpStat(psr, freq)
        
        # check to see if larger than in real data
        if fpstat > fstat_ref:
            count += 1

    # now get detection probability
    detProb = count/nreal

    print freq, h, detProb

    return detProb - 0.95

def marginalizedPulsarPhaseLike(psr, theta, phi, phase, inc, psi, freq, h, maximize=False):
    """ 
    Compute the log-likelihood marginalized over pulsar phases

    @param psr: List of pulsar object instances
    @param theta: GW polar angle [radian]
    @param phi: GW azimuthal angle [radian]
    @param phase: Initial GW phase [radian]
    @param inc: GW inclination angle [radian]
    @param psi: GW polarization angle [radian]
    @param freq: GW initial frequency [Hz]
    @param h: GW strain
    @param maximize: Option to maximize over pulsar phases instead of marginalize

    """

    # get number of pulsars
    npsr = len(psr)
       
    # get c and d
    c = np.cos(phase)
    d = np.sin(phase)

    # construct xi = M**5/3/D and omega
    xi = 0.25 * np.sqrt(5/2) * (np.pi*freq)**(-2/3) * h
    omega = np.pi*freq

    lnlike = 0
    for ct, pp in enumerate(psr):

        # compute relevant inner products
        cip = np.dot(np.cos(2*omega*pp.toas), np.dot(pp.invCov, pp.res)) 
        sip = np.dot(np.sin(2*omega*pp.toas), np.dot(pp.invCov, pp.res))
        N = np.dot(np.cos(2*omega*pp.toas), np.dot(pp.invCov, np.cos(2*omega*pp.toas)))

        # compute fplus and fcross
        fplus, fcross, cosMu = PALutils.createAntennaPatternFuncs(pp, theta, phi)

        # mind you p's and q's
        p = (1+np.cos(inc)**2) * (fplus*np.cos(2*psi) + fcross*np.sin(2*psi))
        q = 2*np.cos(inc) * (fplus*np.sin(2*psi) - fcross*np.cos(2*psi))

        # construct X Y and Z
        X = -xi/omega**(1/3) * (p*sip + q*cip - 0.5*xi/omega**(1/3)*N*c*(p**2+q**2))
        Y = -xi/omega**(1/3) * (q*sip - p*cip - 0.5*xi/omega**(1/3)*N*d*(p**2+q**2))
        Z = xi/omega**(1/3) * ((p*c+q*d)*sip - (p*d-q*c)*cip \
                        -0.5*xi/omega**(1/3)*N*(p**2+q**2))

        # add to log-likelihood
        #print X, Y
        if maximize:
            lnlike += Z + np.sqrt(X**2 + Y**2)
        else:
            lnlike += Z + np.log(ss.iv(0, np.sqrt(X**2 + Y**2)))

    return lnlike

def crossPower(psr, gam=13/3):
    """

    Compute the cross power as defined in Eq 9 and uncertainty of Eq 10 in 
    Demorest et al (2012).

    @param psr: List of pulsar object instances
    @param gam: Power spectral index of GWB

    @return: vector of cross power for each pulsar pair
    @return: vector of cross power uncertainties for each pulsar pair

    """

    # initialization
    npsr = len(psr) 

    # now compute cross power
    rho = []
    sig = []
    xi = []
    for ll in range(npsr):
        for kk in range(ll+1, npsr):
            
            # matrix of time lags
            tm = PALutils.createTimeLags(psr[ll].toas, psr[kk].toas)

            # create cross covariance matrix without overall amplitude A^2
            SIJ = PALutils.createRedNoiseCovarianceMatrix(tm, 1, gam)
            
            # construct numerator and denominator of optimal statistic
            bot = np.trace(np.dot(psr[ll].invCov, np.dot(SIJ, np.dot(psr[kk].invCov, SIJ.T))))
            top = np.dot(psr[ll].res, np.dot(psr[ll].invCov, np.dot(SIJ, \
                        np.dot(psr[kk].invCov, psr[kk].res))))

            # cross correlation and uncertainty
            rho.append(top/bot)
            sig.append(1/np.sqrt(bot))


    return np.array(rho), np.array(sig)

def bigPulsarPhaseJump(x, iter, beta):

    # get old parameters
    q = x.copy()

    # pick pulsar index at random
    ind = np.random.randint(0, npsr, npsr)
    ind = np.unique(ind)
    
    # get relevant parameters
    freq = 10**x[2]
    pdist = x[8+ind]
    pdistErr = np.array([psr[ii].distErr for ii in list(ind)])
    phi = x[1]
    theta = x[0]

    # put pulsar distance in correct units
    pdist *= 1.0267e11  
    pdistErr *= 1.0267e11  

    # get cosMu
    cosMu = np.zeros(len(ind))
    for ii in range(len(ind)):
        tmp1, temp2, cosMu[ii] = PALutils.createAntennaPatternFuncs(psr[ind[ii]], theta, phi)
    
    # construct pulsar phase
    phase_old = 2*np.pi*freq*pdist*(1-cosMu)

    # gaussian jump 
    phase_jump = np.random.randn(np.size(pdist))*pdistErr*freq*(1-cosMu)

    # make jump multiple of 2 pi
    phase_jump = np.array([int(phase_jump[ii]) \
                    for ii in range(np.size(pdist))])

    # new phase
    phase_new = phase_old + 2*np.pi*phase_jump

    # solve for new pulsar distances from phase_new
    L_new = phase_new/(2*np.pi*freq*(1-cosMu))

    # convert back to Kpc
    L_new /= 1.0267e11  

    q[8+ind] = L_new

    return q

def smallPulsarPhaseJump(x, iter, beta):

    # get old parameters
    q = x.copy()

    # jump size
    sigma = np.sqrt(beta) * 0.5

    # pick pulsar index at random
    ind = np.random.randint(0, npsr, npsr)
    ind = np.unique(ind)
    
    # get relevant parameters
    freq = 10**x[2]
    pdist = x[8+ind]
    phi = x[1]
    theta = x[0]

    # put pulsar distance in correct units
    pdist *= 1.0267e11  

    # get cosMu
    cosMu = np.zeros(len(ind))
    for ii in range(len(ind)):
        tmp1, temp2, cosMu[ii] = PALutils.createAntennaPatternFuncs(psr[ind[ii]], theta, phi)


    # construct pulsar phase
    phase_old = 2*np.pi*freq*pdist*(1-cosMu)

    # gaussian jump 
    phase_new = phase_old + np.random.randn(np.size(pdist))*sigma

    # solve for new pulsar distances from phase_new
    L_new = phase_new/(2*np.pi*freq*(1-cosMu))

    # convert back to Kpc
    L_new /= 1.0267e11 

    q[8+ind] = L_new

    return q

def loglike(x):

    tstart = time.time()

    theta = np.arccos(x[0])
    phi = x[1]
    f = 10**x[2]
    h = 10**x[3]
    psi = x[4]
    inc = np.arccos(x[5])
    phase = x[6]

    # get pulsar phase
    pphase = np.zeros(npsr)
    for ii in range(npsr):
        pphase[ii] = x[(ndim-npsr) + ii]

    # pick a distance and mass from the strain. Doesnt matter since non-evolving
    mc = 5e8
    dist = 4 * np.sqrt(2/5) * (mc*4.9e-6)**(5/3) * (np.pi*f)**(2/3) / h
    dist /= 1.0267e14

    loglike = 0
    for ct, p in enumerate(psr):
    

        # make waveform with no frequency evolution
        s = PALutils.createResiduals(p, theta, phi, mc, dist, f, phase, psi, inc,\
                                     pphase=pphase[ct], evolve=False)

        diff = p.res - s
        loglike += -0.5 * logdetTerm[ct]
        loglike += -0.5 * np.dot(diff, np.dot(p.invCov, diff))

    #print 'Evaluation time = {0} s'.format(time.time() - tstart)

    if np.isnan(loglike):
        print 'NaN log-likelihood. Not good...'
        return -np.inf
    else:
        return loglike

def HTFSingleInd(
    psr,
    F,
    proj,
    SS,
    rho,
    efac,
    equad,
    gwtheta,
    gwphi,
    mc,
    dist,
    fgw,
    phase0,
    psi,
    inc,
    pphase=None,
    pdist=None,
    evolve=True,
    psrTerm=True,
    phase_approx=False,
):
    """
    Lentati marginalized likelihood function only including efac and equad
    and power law coefficients

    @param psr: Pulsar class
    @param F: Fourier design matrix constructed in PALutils
    @param proj: Projection operator from white noise
    @param SS: Diagonalized white noise matrix
    @param rho: Power spectrum coefficients
    @param efac: constant multipier on error bar covaraince matrix term
    @param equad: Additional white noise added in quadrature to efac
    @param gwtheta: Polar angle of GW source in celestial coords [radians]
    @param gwphi: Azimuthal angle of GW source in celestial coords [radians]
    @param mc: Chirp mass of SMBMB [solar masses]
    @param dist: Luminosity distance to SMBMB [Mpc]
    @param fgw: Frequency of GW (twice the orbital frequency) [Hz]
    @param phase0: Initial Phase of GW source [radians]
    @param psi: Polarization of GW source [radians]
    @param inc: Inclination of GW source [radians]
    @param pdist: Pulsar distance to use other than those in psr [kpc]
    @param pphase: Use pulsar phase to determine distance [radian]
    @param psrTerm: Option to include pulsar term [boolean] 
    @param evolve: Option to exclude evolution [boolean]



    @return: LogLike: loglikelihood

    """
    # make waveform with no frequency evolution
    s = PALutils.createResiduals(
        psr,
        gwtheta,
        gwphi,
        mc,
        dist,
        fgw,
        phase0,
        psi,
        inc,
        pphase=pphase,
        evolve=evolve,
        pdist=pdist,
        psrTerm=psrTerm,
        phase_approx=phase_approx,
    )

    diff = np.dot(proj, (psr.res - s))

    # compute total time span of data
    Tspan = psr.toas.max() - psr.toas.min()

    # compute d
    d = np.dot(F.T, diff / (efac * SS + equad ** 2))

    # compute Sigma
    N = 1 / (efac * SS + equad ** 2)
    right = (N * F.T).T
    FNF = np.dot(F.T, right)

    arr = np.zeros(2 * len(rho))
    ct = 0
    for ii in range(0, 2 * len(rho), 2):
        arr[ii] = rho[ct]
        arr[ii + 1] = rho[ct]
        ct += 1

    Sigma = FNF + np.diag(1 / arr)

    # cholesky decomp for second term in exponential
    cf = sl.cho_factor(Sigma)
    expval2 = sl.cho_solve(cf, d)
    logdet_Sigma = np.sum(2 * np.log(np.diag(cf[0])))

    dtNdt = np.sum(diff ** 2 / (efac * SS + equad ** 2))

    logdet_Phi = np.sum(np.log(arr))

#.........这里部分代码省略.........

# make sure all pulsar have same reference time
tt = []
for p in psr:
    tt.append(np.min(p.toas))

# find reference time
tref = np.min(tt)

# now scale pulsar time
for p in psr:
    p.toas -= tref

# get G matrices
for p in psr:
    p.G = PALutils.createGmatrix(p.dmatrix)

# run Fp statistic to determine starting frequency
if args.null == False:
    print "Running initial Fpstat search"
    fsearch = np.logspace(-9, -7, 1000)
    fpstat = np.zeros(len(fsearch))
    for ii in range(len(fsearch)):
        fpstat[ii] = PALLikelihoods.fpStat(psr, fsearch[ii])

    # determine maximum likelihood frequency
    fmaxlike = fsearch[np.argmax(fpstat)]
    print "Maximum likelihood from f-stat search = {0}\n".format(fmaxlike)

# get Tmax
Tmax = np.max([p.toas.max() - p.toas.min() for p in psr])

    def __init__(self,pulsargroup, addNoise=False, addGmatrix=True):


        # loop though keys in pulsargroup and fill in psr attributes that are needed for GW analysis
        self.dist = None
        self.distErr = None
        self.fH = None

        for key in pulsargroup:

            # look for TOAs
            if key == "TOAs":
                self.toas = pulsargroup[key].value

            # residuals
            elif key == "residuals":
                self.res = pulsargroup[key].value

            # toa error bars
            elif key == "toaErr":
                self.err = pulsargroup[key].value
            
            # frequencies in Hz
            elif key == "freqs":
                self.freqs = pulsargroup[key].value
            
            # design matrix
            elif key == "designmatrix":
                self.dmatrix = pulsargroup[key].value
                self.ntoa, self.nfit = self.dmatrix.shape
                if addGmatrix:
                    self.G = PALutils.createGmatrix(self.dmatrix)
            
            # design matrix
            elif key == "pname":
                self.name = pulsargroup[key].value
            
            # pulsar distance in kpc
            elif key == "dist":
                self.dist = pulsargroup[key].value 
            
            # pulsar distance uncertainty in kpc
            elif key == "distErr":
                self.distErr = pulsargroup[key].value 

            # right ascension and declination
            elif key == 'tmp_name':
                par_names = list(pulsargroup[key].value)
                for ct,name in enumerate(par_names):

                    # right ascension and phi
                    if name == "RAJ":
                        self.ra = pulsargroup["tmp_valpost"].value[ct]
                        self.phi = self.ra
                    
                    # right ascension
                    if name == "DECJ":
                        self.dec = pulsargroup["tmp_valpost"].value[ct]
                        self.theta = np.pi/2 - self.dec

            # inverse covariance matrix
            elif key == "invCov":
                if addNoise:
                    self.invCov = pulsargroup[key].value


            ## noise parameters ##

            elif key == "Amp":

                self.Amp = pulsargroup[key].value 
            
            # red noise spectral
            elif key == "gam":
                self.gam = pulsargroup[key].value 
            
            # efac
            elif key == "efac":
                self.efac = pulsargroup[key].value 
            
            # equad
            elif key == "equad":
                self.equad = pulsargroup[key].value 
            
            # fH
            elif key == "fH":
                self.fH = pulsargroup[key].value 

            # pulsar distance uncertainty in kpc
            elif key == "distErr":
                self.distErr = pulsargroup[key].value 

        if self.dist is None:
            print 'WARNING: No distance info, using d = 1 kpc'
            self.dist = 1.0

        if self.distErr is None:
            print 'WARNING: No distance error info, using sigma_d = 0.1 kpc'
            self.distErr = 0.1

    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)
#.........这里部分代码省略.........