Python PALutils.ptSum方法代码示例

# 需要导入模块: import PALutils [as 别名]
# 或者: from PALutils import ptSum [as 别名]
# carry out Fp search
if args.fpFlag:

    print 'Beginning Fp Search with {0} pulsars, with frequency range {1} -- {2}'.format(npsr, f[0], f[-1])
    
    fpstat = np.zeros(args.nfreqs)
    for ii in range(args.nfreqs):

        fpstat[ii] = PALLikelihoods.fpStat(psr, f[ii])


    print 'Done Search. Computing False Alarm Probability'
    
    # single template FAP
    pf = np.array([PALutils.ptSum(npsr, fpstat[ii]) for ii in range(np.alen(f))])

    # get total false alarm probability with trials factor
    pfT = 1 - (1-pf)**np.alen(f)

    # write results to file
    if not os.path.exists(args.outDir):
        os.makedirs(args.outDir)

    # get filename from hdf5 file
    fname = args.outDir + '/' + args.h5file.split('/')[-1].split('.')[0] + '.txt'
    fout = open(fname, 'w')
    print 'Writing results to file {0}'.format(fname)
    for ii in range(np.alen(f)):

        fout.write('%g %g %g\n'%(f[ii], fpstat[ii], pfT[ii]))

# 需要导入模块: import PALutils [as 别名]
# 或者: from PALutils import ptSum [as 别名]
def upperLimitFunc(h, fstat_ref, freq, nreal, theta=None, phi=None, detect=False, \
                  dist=None):
    """
    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(-1, 1))
        #gwpsi = np.random.uniform(-np.pi/4, np.pi/4)
        gwpsi = np.random.uniform(0, np.pi)

        # 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
        
        # check for fixed sky location
        if theta is not None:
            gwtheta = theta
        if phi is not None:
            gwphi = phi
        if dist is not None:
            gwdist = dist
            gwmc = ((gwdist*1.0267e14)/4/np.sqrt(2/5)/(np.pi*freq)**(2/3)*h)**(3/5)/4.9e-6
        
        
        # 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
            noise = np.dot(L[ct], np.random.randn(L[ct].shape[0]))
            p.res = np.dot(R[ct], noise+inducedRes)

        # compute f-statistic
        fpstat = PALLikelihoods.fpStat(psr, freq)
        
        # check to see if larger than in real data
        if detect:
            if PALutils.ptSum(npsr, fpstat) < 1e-4:
                count += 1
        else:
            if fpstat > fstat_ref:
                count += 1

    # now get detection probability
    detProb = count/nreal
    
    if args.dist:
        print '%e %e %f\n'%(freq, gwmc, detProb)
    else:
        print freq, h, detProb

    return detProb - 0.95