Python PALutils.createAntennaPatternFuncs方法代码示例

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

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

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

# 需要导入模块: import PALutils [as 别名]
# 或者: from PALutils import createAntennaPatternFuncs [as 别名]
def modelIndependentFullPTASinglSource(psr, proj, s, f, theta, phi, rho, kappa, efac, equad, ORF):
    """
    Model Independent single source testing function

    """
    tstart = time.time()
    
    # get the number of modes, should be the same for all pulsars
    nmode = len(rho)
    npsr = len(psr)

    # get F matrices for all pulsars at given frequency
    F = [np.array([np.sin(2*np.pi*f*p.toas), np.cos(2*np.pi*f*p.toas)]).T for p in psr]

    F = [np.dot(proj[ii], F[ii]) for ii in range(len(proj))]

    loglike1 = 0
    FtNF = []
    for ct,p in enumerate(psr):
    
        # compute d
        if ct == 0:
            d = np.dot(F[ct].T, p.res/(efac[ct]*s[ct] + equad[ct]**2))
        else:
            d = np.append(d, np.dot(F[ct].T, p.res/(efac[ct]*s[ct] + equad[ct]**2)))

        # compute FT N F
        N = 1/(efac[ct]*s[ct] + equad[ct]**2)
        right = (N*F[ct].T).T
        FtNF.append(np.dot(F[ct].T, right))
        
        # log determinant of N
        logdet_N = np.sum(np.log(efac[ct]*s[ct] + equad[ct]**2))

        # triple produce in likelihood function
        dtNdt = np.sum(p.res**2/(efac[ct]*s[ct] + equad[ct]**2))

        loglike1 += -0.5 * (logdet_N + dtNdt)

    # construct elements of sigma array
    sigdiag = []
    sigoffdiag = []
    fplus = np.zeros(npsr)
    fcross = np.zeros(npsr)
    for ii in range(npsr):
        fplus[ii], fcross[ii], cosMu = PALutils.createAntennaPatternFuncs(psr[ii], theta, phi)
        tot = np.zeros(2*nmode)
        offdiag = np.zeros(2*nmode)

        # off diagonal terms
        offdiag[0::2] = 10**rho 
        offdiag[1::2] = 10**rho

        # diagonal terms
        tot[0::2] = 10**rho
        tot[1::2] = 10**rho

        # add in individual red noise
        if len(kappa[ii]) > 0:
            tot[0::2][0:len(kappa[ii])] += 10**kappa[ii]
            tot[1::2][0:len(kappa[ii])] += 10**kappa[ii]
        
        # fill in lists of arrays
        sigdiag.append(tot)
        sigoffdiag.append(offdiag)

    tstart2 = time.time()

    # compute Phi inverse from Lindley's code
    smallMatrix = np.zeros((2*nmode, npsr, npsr))
    for ii in range(npsr):
        for jj in range(ii,npsr):

            if ii == jj:
                smallMatrix[:,ii,jj] = ORF[ii,jj] * sigdiag[jj] * (fplus[ii]**2 + fcross[ii]**2)
            else:
                smallMatrix[:,ii,jj] = ORF[ii,jj] * sigoffdiag[jj] * (fplus[ii]*fplus[jj] + fcross[ii]*fcross[jj])
                smallMatrix[:,jj,ii] = smallMatrix[:,ii,jj]


    # invert them
    logdet_Phi = 0
    for ii in range(2*nmode):
        L = sl.cho_factor(smallMatrix[ii,:,:])
        smallMatrix[ii,:,:] = sl.cho_solve(L, np.eye(npsr))
        logdet_Phi += np.sum(2*np.log(np.diag(L[0])))

    # now fill in real covariance matrix
    Phi = np.zeros((2*npsr*nmode, 2*npsr*nmode))
    for ii in range(npsr):
        for jj in range(ii,npsr):
            for kk in range(0,2*nmode):
                Phi[kk+ii*2*nmode,kk+jj*2*nmode] = smallMatrix[kk,ii,jj]
    
    # symmeterize Phi
    Phi = Phi + Phi.T - np.diag(np.diag(Phi))
            
    # compute sigma
    Sigma = sl.block_diag(*FtNF) + Phi

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

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

    """

    tstart = time.time()

    # get number of pulsars
    npsr = len(psr)
       
    # 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
    
    # get a values from Ellis et al 2012
    a1 = xi * ((1+np.cos(inc)**2)*np.cos(phase)*np.cos(2*psi) + \
               2*np.cos(inc)*np.sin(phase)*np.sin(2*psi))
    a2 = -xi * ((1+np.cos(inc)**2)*np.sin(phase)*np.cos(2*psi) - \
                2*np.cos(inc)*np.cos(phase)*np.sin(2*psi))
    a3 = xi * ((1+np.cos(inc)**2)*np.cos(phase)*np.sin(2*psi) - \
               2*np.cos(inc)*np.sin(phase)*np.cos(2*psi))
    a4 = -xi * ((1+np.cos(inc)**2)*np.sin(phase)*np.sin(2*psi) + \
                2*np.cos(inc)*np.cos(phase)*np.cos(2*psi))

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

        tstartip = time.time()

        # compute relevant inner products
        N1 = np.dot(np.cos(2*omega*pp.toas), np.dot(pp.invCov, pp.res)) 
        N2 = np.dot(np.sin(2*omega*pp.toas), np.dot(pp.invCov, pp.res))
        M11 = np.dot(np.sin(2*omega*pp.toas), np.dot(pp.invCov, np.sin(2*omega*pp.toas)))
        M22 = np.dot(np.cos(2*omega*pp.toas), np.dot(pp.invCov, np.cos(2*omega*pp.toas)))
        M12 = np.dot(np.cos(2*omega*pp.toas), np.dot(pp.invCov, np.sin(2*omega*pp.toas)))

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

        # mind your p's and q's
        p = fplus*a1 + fcross*a3
        q = fplus*a2 + fcross*a4

        # constuct multipliers of pulsar phase terms
        X = p*N1 + q*N2 + p**2*M11 + q**2*M22 + 2*p*q*M12
        Y = p*N1 + q*N2 + 2*p**2*M11 + 2*q**2*M22 + 4*p*q*M12
        Z = p*N2 - q*N1 + 2*(p**2-q**2)*M12 - 2*p*q*(M11-M22)
        W = q**2*M11 + p**2*M22 -2*p*q*M12
        V = p*q*(M11-M22) - (p**2-q**2)*M12
        
        #print X, Y, Z, W, V
        tip += (time.time() - tstartip)

        tstartint = time.time()

        # find the maximum of argument of exponential function
        phip = np.linspace(0, 2*np.pi, 10000)
        arg = X - Y*np.cos(phip) + Z*np.sin(phip) + W*np.sin(phip)**2 + 2*V*np.cos(phip)*np.sin(phip)
        maxarg = np.max(arg)

        if maximize:
            tmax += maxarg

        else:

            # define integrand for numerical integration
            f = lambda phi: np.exp(X - Y*np.cos(phi) + Z*np.sin(phi) + \
                    W*np.sin(phi)**2 + 2*V*np.cos(phi)*np.sin(phi) - maxarg)

            # do numerical integration
            integral = integrate.quad(f, 0, 2*np.pi)[0]
            lnlike += maxarg + np.log(integral)

            tint += (time.time() - tstartint)

    print 'Loglike = {0}'.format(lnlike)
    print 'Total Evaluation Time = {0} s'.format(time.time() - tstart)
    print 'Total inner product evaluation Time = {0} s'.format(tip)
    print 'Total Integration Time = {0} s\n'.format(tint)

    if maximize:
        lnlike = tmax

    return lnlike

# 需要导入模块: import PALutils [as 别名]
# 或者: from PALutils import createAntennaPatternFuncs [as 别名]
def covarianceJumpProposal(x, iter, beta):

    # get scale
    scale = 1

    # medium size jump every 1000 steps
    if np.random.rand() < 1/1000:
        scale = 5
    
    # large size jump every 10000 steps
    if np.random.rand() < 1/10000:
        scale = 50

    # get parmeters in new diagonalized basis
    y = np.dot(U.T, x)

    # make correlated componentwise adaptive jump
    if iter < 100000:
        ind = np.unique(np.random.randint(0, 8, 1))
    else:
        ind = np.unique(np.random.randint(0, 8, ndim))

    neff = len(ind)
    cd = 2.4 * np.sqrt(1/beta) / np.sqrt(neff) * scale
    #y[ind] = y[ind] + np.random.randn(neff) * cd * np.sqrt(S[ind])
    q = x.copy()
    q[ind] = q[ind] + np.random.randn(neff) * cd * np.sqrt(np.diag(cov[ind,ind]))
    #q = np.dot(U, y)

    # need to make sure that we keep the pulsar phase constant, plus small offset

    # get parameters before jump
    omega0 = 2*np.pi*10**x[2]
    L0 = x[8:] * 1.0267e11 
    phi0 = x[1]
    theta0 = x[0]
    
    omega1 = 2*np.pi*10**q[2]
    L1 = q[8:] * 1.0267e11 
    phi1 = q[1]
    theta1 = q[0]

    # get cosMu
    cosMu0 = np.zeros(npsr)
    cosMu1 = np.zeros(npsr)
    for ii in range(npsr):
        tmp1, temp2, cosMu0[ii] = PALutils.createAntennaPatternFuncs(psr[ii], theta0, phi0)
        tmp1, temp2, cosMu1[ii] = PALutils.createAntennaPatternFuncs(psr[ii], theta1, phi1)
    
    # construct new pulsar distances to keep the pulsar phases constant
    sigma = np.sqrt(1/beta) * 0.0 * np.random.randn(npsr)
    phase0 = omega0*L0*(1-cosMu0)
    phase1 = omega1*L1*(1-cosMu1)
    L_new = L1 - (np.mod(phase1, 2*np.pi) - np.mod(phase0, 2*np.pi) + sigma)/(omega1*(1-cosMu1)) 

    #phasenew = omega1*L_new*(1-cosMu1)
    #print 'New phase = ', np.mod(phasenew ,2*np.pi) 
    #print 'Old phase = ', np.mod(phase0,2*np.pi)
    #print 'New L = ', L_new 
    #print 'Old L = ', L0, '\n'

    # convert back to Kpc
    L_new /= 1.0267e11 

    q[8:] = L_new

    #print x[0]
    #print q[0], '\n'

    return q