test_pta.py 文件源码

def test_parameterized_orf(self):
        T1 = 3.16e8
        pl = utils.powerlaw(log10_A=parameter.Uniform(-18,-12),
                            gamma=parameter.Uniform(1,7))
        orf = hd_orf_generic(a=parameter.Uniform(0,5),
                             b=parameter.Uniform(0,5),
                             c=parameter.Uniform(0,5))
        rn = gp_signals.FourierBasisGP(spectrum=pl, Tspan=T1, components=30)
        crn = gp_signals.FourierBasisCommonGP(spectrum=pl, orf=orf,
                                              components=30, name='gw',
                                              Tspan=T1)

        model = rn + crn
        pta = model(self.psrs[0]) + model(self.psrs[1])

        lA1, gamma1 = -13, 1e-15
        lA2, gamma2 = -13.3, 1e-15
        lAc, gammac = -13.1, 1e-15
        a, b, c = 1.9, 0.4, 0.23

        params = {'gw_log10_A': lAc, 'gw_gamma': gammac,
                  'gw_a': a, 'gw_b':b, 'gw_c':c,
                  'B1855+09_log10_A': lA1, 'B1855+09_gamma': gamma1,
                  'J1909-3744_log10_A': lA2, 'J1909-3744_gamma': gamma2}

        phi = pta.get_phi(params)
        phiinv = pta.get_phiinv(params)

        F1, f1 = utils.createfourierdesignmatrix_red(
            self.psrs[0].toas, nmodes=30, Tspan=T1)
        F2, f2 = utils.createfourierdesignmatrix_red(
            self.psrs[1].toas, nmodes=30, Tspan=T1)

        msg = 'F matrix incorrect'
        assert np.allclose(pta.get_basis(params)[0], F1, rtol=1e-10), msg
        assert np.allclose(pta.get_basis(params)[1], F2, rtol=1e-10), msg

        nftot = 120
        phidiag = np.zeros(nftot)
        phit = np.zeros((nftot, nftot))

        phidiag[:60] = utils.powerlaw(f1, lA1, gamma1)
        phidiag[:60] += utils.powerlaw(f1, lAc, gammac)
        phidiag[60:] = utils.powerlaw(f2, lA2, gamma2)
        phidiag[60:] += utils.powerlaw(f2, lAc, gammac)

        phit[np.diag_indices(nftot)] = phidiag
        orf = hd_orf_generic(self.psrs[0].pos, self.psrs[1].pos, a=a, b=b, c=c)
        spec = utils.powerlaw(f1, log10_A=lAc, gamma=gammac)
        phit[:60, 60:] = np.diag(orf*spec)
        phit[60:, :60] = phit[:60, 60:]

        msg = '{} {}'.format(np.diag(phi), np.diag(phit))
        assert np.allclose(phi, phit, rtol=1e-15, atol=1e-17), msg
        msg = 'PTA Phi inverse is incorrect {}.'.format(params)
        assert np.allclose(phiinv, np.linalg.inv(phit),
                           rtol=1e-15, atol=1e-17), msg