Python PALutils.createResiduals方法代码示例

# 需要导入模块: import PALutils [as 别名]
# 或者: from PALutils import createResiduals [as 别名]
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

# 需要导入模块: import PALutils [as 别名]
# 或者: from PALutils import createResiduals [as 别名]
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

# 需要导入模块: import PALutils [as 别名]
# 或者: from PALutils import createResiduals [as 别名]
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

# 需要导入模块: import PALutils [as 别名]
# 或者: from PALutils import createResiduals [as 别名]
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))

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

# 需要导入模块: import PALutils [as 别名]
# 或者: from PALutils import createResiduals [as 别名]
def upperLimitFunc(h):
    """
    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)

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

        # 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 and refit for all pulsars
        for ct,p in enumerate(psr):
            inducedRes = PALutils.createResiduals(p, gwtheta, gwphi, gwmc, gwdist, \
                            freq, gwphase, gwpsi, gwinc)
 
            # create simulated data set
            noise = np.dot(L[ct], np.random.randn(L[ct].shape[0]))
            pp[ct].stoas[:] -= pp[ct].residuals()/86400
            pp[ct].stoas[:] += np.longdouble(np.dot(RQ[ct], noise)/86400)
            pp[ct].stoas[:] += np.longdouble(np.dot(RQ[ct], inducedRes)/86400)

            # refit
            pp[ct].fit(iters=3)

            # replace residuals in pulsar object
            p.res = pp[ct].residuals()

            print p.name, p.rms()*1e6

        # 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 h, detProb

    return detProb - 0.95

# 需要导入模块: import PALutils [as 别名]
# 或者: from PALutils import createResiduals [as 别名]
    if args.nopterm:
        print 'Not including pulsar term in single source waveform!'
        pterm = False
    else:
        pterm = True
    
    if args.snr is not None:
          
        print 'Scaling distance to give SNR = {0}'.format(args.snr)

        snr2 = 0
        args.gwdist = 1
        for ct, p in enumerate(pp):

            inducedRes = (PALutils.createResiduals(psr[ct], np.pi/2-args.gwdec, args.gwra, \
                            args.gwchirpmass, args.gwdist, args.gwfreq, args.gwphase, \
                            args.gwpolarization, args.gwinc, psrTerm=pterm))

            # compute snr
            snr2 += PALutils.calculateMatchedFilterSNR(psr[ct], inducedRes, inducedRes)**2

        # get total snr
        snr = np.sqrt(snr2)

        # scale distance appropiately
        args.gwdist = snr/args.snr

        print 'Scaled GW distance = {0} for SNR = {1}'.format(args.gwdist, args.snr)
    
    # make residuals
    for ct, p in enumerate(pp):