From 1a1c9d9955c176d0cfbb1c40e4875eecbd870e0b Mon Sep 17 00:00:00 2001 From: Stephen Taylor Date: Thu, 30 May 2019 16:02:43 -0700 Subject: [PATCH] Editing utils.py to patch py27 functionality (#193) * adding ability to have user-defined modes in CommonGP * fixing lint issues * adding dynamic partial files * adding code for dynamic BayesEphem * fixing ephemeris tests * fixing ephemeris tests * fixing ephemeris tests * fixing to allow py27 support * running travis tests on py2 and py3 * changing numpy requirements * changing python requirements * patching to avoid crashing when python2 is used with gp coefficients * final lint corrections --- .travis.yml | 1 + enterprise/signals/gp_signals.py | 21 +- enterprise/signals/utils.py | 7 +- requirements.txt | 2 +- tests/test_gp_coefficients.py | 406 ++++++++++++++++--------------- 5 files changed, 238 insertions(+), 199 deletions(-) diff --git a/.travis.yml b/.travis.yml index 5d7697e1..38980229 100644 --- a/.travis.yml +++ b/.travis.yml @@ -2,6 +2,7 @@ language: python cache: pip python: + - "2.7.15" - "3.6" before_install: diff --git a/enterprise/signals/gp_signals.py b/enterprise/signals/gp_signals.py index e34538d2..87adff7d 100644 --- a/enterprise/signals/gp_signals.py +++ b/enterprise/signals/gp_signals.py @@ -10,6 +10,8 @@ import math import itertools import functools +import platform +import logging import numpy as np @@ -21,6 +23,12 @@ from enterprise.signals.selections import Selection from enterprise.signals.utils import KernelMatrix +pyv3 = platform.python_version().split('.')[0] == '3' + +logging.basicConfig(format='%(levelname)s: %(name)s: %(message)s', + level=logging.INFO) +logger = logging.getLogger(__name__) + def BasisGP(priorFunction, basisFunction, coefficients=False, combine=True, selection=Selection(selections.no_selection), @@ -39,6 +47,9 @@ def __init__(self, psr): self.name = self.psrname + '_' + self.signal_id self._do_selection(psr, priorFunction, basisFunction, coefficients, selection) + if coefficients and not pyv3: + msg = 'GP coefficients compatible only with Python 3' + logger.warning(msg) def _do_selection(self, psr, priorfn, basisfn, coefficients, selection): @@ -60,7 +71,7 @@ def _do_selection(self, psr, priorfn, basisfn, coefficients, self._bases[key]._params.values()): self._params[par.name] = par - if coefficients: + if coefficients and pyv3: # we can only create GPCoefficients parameters if the basis # can be constructed with default arguments # (and does not change size) @@ -116,7 +127,7 @@ def _construct_basis(self, params={}): # this class does different things (and gets different method # definitions) if the user wants it to model GP coefficients # (e.g., for a hierarchical likelihood) or if they do not - if coefficients: + if coefficients and pyv3: def _get_coefficient_logprior(self, key, c, **params): self._construct_basis(params) @@ -230,7 +241,7 @@ class TimingModel(BaseClass): signal_name = 'linear timing model' signal_id = name + '_svd' if use_svd else name - if coefficients: + if coefficients and pyv3: def _get_coefficient_logprior(self, key, c, **params): # MV: probably better to avoid this altogether # than to use 1e40 as in get_phi @@ -296,7 +307,7 @@ def __init__(self, psr): self._psrpos = psr.pos - if coefficients: + if coefficients and pyv3: self._construct_basis() chain = itertools.chain(self._prior._params.values(), @@ -326,7 +337,7 @@ def basis_params(self): def _construct_basis(self, params={}): self._basis, self._labels = self._bases(params=params) - if coefficients: + if coefficients and pyv3: def _get_coefficient_logprior(self, c, **params): # MV: for correlated GPs, the prior needs to use # the coefficients for all GPs together; diff --git a/enterprise/signals/utils.py b/enterprise/signals/utils.py index 9d881a8c..c6446a00 100644 --- a/enterprise/signals/utils.py +++ b/enterprise/signals/utils.py @@ -1178,8 +1178,13 @@ def createfourierdesignmatrix_physicalephem(toas, planetssb, pos_t, c = np.zeros(6) c[i] = dpar + #Fl.append(physical_ephem_delay(toas, planetssb, pos_t, + # **{parname: c}, **oa)/dpar) + kwarg_dict = {parname: c} + kwarg_dict.update(oa) Fl.append(physical_ephem_delay(toas, planetssb, pos_t, - **{parname: c}, **oa)/dpar) + **kwarg_dict)/dpar) + Phil.append(ppar) return np.array(Fl).T.copy(), np.array(Phil) diff --git a/requirements.txt b/requirements.txt index f4ffe34d..8cb558e7 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,4 @@ -numpy==1.16.0 +numpy==1.16.3 scipy==1.2.0 ephem>=3.7.6.0 Cython==0.28.5 diff --git a/tests/test_gp_coefficients.py b/tests/test_gp_coefficients.py index 5f4a42f1..04fb7a16 100644 --- a/tests/test_gp_coefficients.py +++ b/tests/test_gp_coefficients.py @@ -12,6 +12,8 @@ import math import unittest import numpy as np +import platform +import logging from tests.enterprise_test_data import datadir from enterprise.pulsar import Pulsar @@ -24,6 +26,12 @@ from enterprise.signals import signal_base from enterprise.signals import utils +pyv3 = platform.python_version().split('.')[0] == '3' + +logging.basicConfig(format='%(levelname)s: %(name)s: %(message)s', + level=logging.INFO) +logger = logging.getLogger(__name__) + @signal_base.function def create_quant_matrix(toas, dt=1): @@ -52,226 +60,240 @@ def setUpClass(cls): cls.psr2 = Pulsar(datadir + '/B1937+21_NANOGrav_9yv1.gls.par', datadir + '/B1937+21_NANOGrav_9yv1.tim') - def test_ephemeris(self): - """Test physical-ephemeris delay, made three ways: from marginalized - GP, from coefficient-based GP, from deterministic model.""" - - ef = white_signals.MeasurementNoise(efac=parameter.Uniform(0.1, 5.0)) - - eph = gp_signals.FourierBasisCommonGP_physicalephem( - sat_orb_elements=None) - - ephc = gp_signals.FourierBasisCommonGP_physicalephem( - sat_orb_elements=None, coefficients=True) - - ephd = deterministic_signals.PhysicalEphemerisSignal( - inc_saturn_orb=False) - - model = ef + eph - modelc = ef + ephc - modeld = ef + ephd + if not pyv3: + msg = 'GP coefficients compatible only with Python 3. Skipping tests' + logger.warning(msg) + pass + else: + def test_ephemeris(self): + """Test physical-ephemeris delay, made three ways: from + marginalized GP, from coefficient-based GP, from + deterministic model.""" - pta = signal_base.PTA([model(self.psr), model(self.psr2)]) - ptac = signal_base.PTA([modelc(self.psr), modelc(self.psr2)]) - ptad = signal_base.PTA([modeld(self.psr), modeld(self.psr2)]) + ef = white_signals.MeasurementNoise( + efac=parameter.Uniform(0.1, 5.0)) - cf = 1e-3 * np.random.randn(11) - cf[0] = 1e-5 # this is more sensitive to linearity + eph = gp_signals.FourierBasisCommonGP_physicalephem( + sat_orb_elements=None) - bs = pta.get_basis() - da = [np.dot(bs[0], cf), np.dot(bs[1], cf)] + ephc = gp_signals.FourierBasisCommonGP_physicalephem( + sat_orb_elements=None, coefficients=True) - params = {'B1855+09_efac': 1, 'B1937+21_efac': 1, - 'B1855+09_physicalephem_gp_coefficients': cf, - 'B1937+21_physicalephem_gp_coefficients': cf} - db = ptac.get_delay(params=params) + ephd = deterministic_signals.PhysicalEphemerisSignal( + inc_saturn_orb=False) - dparams = {'B1855+09_efac': 1, 'B1937+21_efac': 1, - 'frame_drift_rate': cf[0], - 'd_jupiter_mass': cf[1], 'd_saturn_mass': cf[2], - 'd_uranus_mass': cf[3], 'd_neptune_mass': cf[4], - 'jup_orb_elements': cf[5:]} - dc = ptad.get_delay(params=dparams) + model = ef + eph + modelc = ef + ephc + modeld = ef + ephd - msg = "Reconstructed ephemeris signals differ!" + pta = signal_base.PTA([model(self.psr), model(self.psr2)]) + ptac = signal_base.PTA([modelc(self.psr), modelc(self.psr2)]) + ptad = signal_base.PTA([modeld(self.psr), modeld(self.psr2)]) - assert np.allclose(da[0], db[0]), msg - assert np.allclose(da[1], db[1]), msg + cf = 1e-3 * np.random.randn(11) + cf[0] = 1e-5 # this is more sensitive to linearity - # we don't expect an exact match since we are linearizing - assert np.allclose(da[0], dc[0], atol=1e-3), msg - assert np.allclose(da[1], dc[1], atol=1e-3), msg + bs = pta.get_basis() + da = [np.dot(bs[0], cf), np.dot(bs[1], cf)] - def test_common_red_noise(self): - """Test of a coefficient-based common GP.""" - pl = utils.powerlaw(log10_A=parameter.Uniform(-18,-12), - gamma=parameter.Uniform(1,7)) + params = {'B1855+09_efac': 1, 'B1937+21_efac': 1, + 'B1855+09_physicalephem_gp_coefficients': cf, + 'B1937+21_physicalephem_gp_coefficients': cf} + db = ptac.get_delay(params=params) - ef = white_signals.MeasurementNoise(efac=parameter.Uniform(0.1, 5.0)) + dparams = {'B1855+09_efac': 1, 'B1937+21_efac': 1, + 'frame_drift_rate': cf[0], + 'd_jupiter_mass': cf[1], 'd_saturn_mass': cf[2], + 'd_uranus_mass': cf[3], 'd_neptune_mass': cf[4], + 'jup_orb_elements': cf[5:]} + dc = ptad.get_delay(params=dparams) - Tspan = (max(self.psr.toas.max(), self.psr2.toas.max()) - - min(self.psr.toas.max(), self.psr2.toas.max())) + msg = "Reconstructed ephemeris signals differ!" - pl = utils.powerlaw(log10_A=parameter.Uniform(-18,-12), - gamma=parameter.Uniform(1,7)) + assert np.allclose(da[0], db[0]), msg + assert np.allclose(da[1], db[1]), msg - rn = gp_signals.FourierBasisCommonGP(spectrum=pl, orf=utils.hd_orf(), - components=20, Tspan=Tspan) + # we don't expect an exact match since we are linearizing + assert np.allclose(da[0], dc[0], atol=1e-3), msg + assert np.allclose(da[1], dc[1], atol=1e-3), msg - model = ef + rn + def test_common_red_noise(self): + """Test of a coefficient-based common GP.""" + pl = utils.powerlaw(log10_A=parameter.Uniform(-18,-12), + gamma=parameter.Uniform(1,7)) - rnc = gp_signals.FourierBasisCommonGP(spectrum=pl, orf=utils.hd_orf(), - components=20, Tspan=Tspan, - coefficients=True) + ef = white_signals.MeasurementNoise( + efac=parameter.Uniform(0.1, 5.0)) - modelc = ef + rnc + Tspan = (max(self.psr.toas.max(), self.psr2.toas.max()) - + min(self.psr.toas.max(), self.psr2.toas.max())) - pta = signal_base.PTA([model(self.psr), model(self.psr2)]) - ptac = signal_base.PTA([modelc(self.psr), modelc(self.psr2)]) + pl = utils.powerlaw(log10_A=parameter.Uniform(-18,-12), + gamma=parameter.Uniform(1,7)) - params = {'B1855+09_efac': 1.0, - 'B1937+21_efac': 1.0, - 'common_fourier_gamma': 5, - 'common_fourier_log10_A': -15} + rn = gp_signals.FourierBasisCommonGP(spectrum=pl, + orf=utils.hd_orf(), + components=20, Tspan=Tspan) - # get GP delays in two different ways + model = ef + rn - cf, cf2 = np.random.randn(40), np.random.randn(40) + rnc = gp_signals.FourierBasisCommonGP(spectrum=pl, + orf=utils.hd_orf(), + components=20, Tspan=Tspan, + coefficients=True) - bs = pta.get_basis(params) - da = [np.dot(bs[0], cf), np.dot(bs[1], cf2)] + modelc = ef + rnc - params.update({'B1855+09_common_fourier_coefficients': cf, - 'B1937+21_common_fourier_coefficients': cf2}) + pta = signal_base.PTA([model(self.psr), model(self.psr2)]) + ptac = signal_base.PTA([modelc(self.psr), modelc(self.psr2)]) - db = ptac.get_delay(params) + params = {'B1855+09_efac': 1.0, + 'B1937+21_efac': 1.0, + 'common_fourier_gamma': 5, + 'common_fourier_log10_A': -15} - msg = 'Implicit and explicit GP delays are different.' - assert np.allclose(da[0], db[0]), msg - assert np.allclose(da[1], db[1]), msg + # get GP delays in two different ways - cpar = [p for p in ptac.params if 'coefficients' in p.name] + cf, cf2 = np.random.randn(40), np.random.randn(40) - def shouldfail(): - return cpar[0].get_logpdf(params) + bs = pta.get_basis(params) + da = [np.dot(bs[0], cf), np.dot(bs[1], cf2)] - self.assertRaises(NotImplementedError, shouldfail) + params.update({'B1855+09_common_fourier_coefficients': cf, + 'B1937+21_common_fourier_coefficients': cf2}) - def test_fourier_red_noise(self): - """Test that implicit and explicit GP delays are the same.""" - # set up signal parameter - pl = utils.powerlaw(log10_A=parameter.Uniform(-18,-12), - gamma=parameter.Uniform(1,7)) - rn = gp_signals.FourierBasisGP(spectrum=pl, components=20) - rnm = rn(self.psr) - - rnc = gp_signals.FourierBasisGP(spectrum=pl, components=20, - coefficients=True) - rnmc = rnc(self.psr) - - # parameters - log10_A, gamma = -14.5, 4.33 - params = {'B1855+09_red_noise_log10_A': log10_A, - 'B1855+09_red_noise_gamma': gamma} - - # get the GP delays in two different ways - cf = np.random.randn(40) - d1 = np.dot(rnm.get_basis(params),cf) - - params.update({'B1855+09_red_noise_coefficients': cf}) - d2 = rnmc.get_delay(params) - - msg = 'Implicit and explicit GP delays are different.' - assert np.allclose(d1, d2), msg - - # np.array cast is needed because we get a KernelArray - phimat = np.array(rnm.get_phi(params)) - pr1 = (-0.5 * np.sum(cf*cf/phimat) - - 0.5 * np.sum(np.log(phimat)) - - 0.5 * len(phimat) * np.log(2*math.pi)) - - cpar = [p for p in rnmc.params if 'coefficients' in p.name][0] - pr2 = cpar.get_logpdf(params=params) - - msg = 'Implicit and explicit GP priors are different.' - assert np.allclose(pr1, pr2), msg - - def test_ecorr_backend(self): - """Test that ecorr-backend signal returns correct values.""" - # set up signal parameter - ecorr = parameter.Uniform(-10, -5) - selection = Selection(selections.by_backend) - ec = gp_signals.EcorrBasisModel(log10_ecorr=ecorr, - selection=selection) - ecm = ec(self.psr) - - ecc = gp_signals.EcorrBasisModel(log10_ecorr=ecorr, - selection=selection, - coefficients=True) - eccm = ecc(self.psr) - - # parameters - ecorrs = [-6.1, -6.2, -6.3, -6.4] - params = {'B1855+09_basis_ecorr_430_ASP_log10_ecorr': ecorrs[0], - 'B1855+09_basis_ecorr_430_PUPPI_log10_ecorr': ecorrs[1], - 'B1855+09_basis_ecorr_L-wide_ASP_log10_ecorr': ecorrs[2], - 'B1855+09_basis_ecorr_L-wide_PUPPI_log10_ecorr': ecorrs[3]} - - fmat = ecm.get_basis(params) - cf = 1e-6 * np.random.randn(fmat.shape[1]) - d1 = np.dot(fmat, cf) - - for key in ecm._keys: - parname = 'B1855+09_basis_ecorr_' + key + '_coefficients' - params[parname] = cf[ecm._slices[key]] - d2 = eccm.get_delay(params) - - msg = 'Implicit and explicit ecorr-basis delays are different.' - assert np.allclose(d1, d2), msg - - def test_formalism(self): - # create marginalized model - ef = white_signals.MeasurementNoise(efac=parameter.Uniform(0.1, 5.0)) - tm = gp_signals.TimingModel() - ec = gp_signals.EcorrBasisModel(log10_ecorr=parameter.Uniform(-10, -5)) - pl = utils.powerlaw(log10_A=parameter.Uniform(-18,-12), - gamma=parameter.Uniform(1,7)) - rn = gp_signals.FourierBasisGP(spectrum=pl, components=10) - model = ef + tm + ec + rn - pta = signal_base.PTA([model(self.psr)]) - - # create hierarchical model - tmc = gp_signals.TimingModel(coefficients=True) - ecc = gp_signals.EcorrBasisModel(log10_ecorr=parameter.Uniform(-10,-5), - coefficients=True) - rnc = gp_signals.FourierBasisGP(spectrum=pl, components=10, - coefficients=True) - modelc = ef + tmc + ecc + rnc - ptac = signal_base.PTA([modelc(self.psr)]) - - ps = {'B1855+09_efac': 1, - 'B1855+09_basis_ecorr_log10_ecorr': -6, - 'B1855+09_red_noise_log10_A': -14, - 'B1855+09_red_noise_gamma': 3} - psc = utils.get_coefficients(pta, ps) - - d1 = ptac.get_delay(psc)[0] - d2 = (np.dot(pta.pulsarmodels[0].signals[1].get_basis(ps), - psc['B1855+09_linear_timing_model_coefficients']) + - np.dot(pta.pulsarmodels[0].signals[2].get_basis(ps), - psc['B1855+09_basis_ecorr_coefficients']) + - np.dot(pta.pulsarmodels[0].signals[3].get_basis(ps), - psc['B1855+09_red_noise_coefficients'])) - - msg = 'Implicit and explicit PTA delays are different.' - assert np.allclose(d1, d2), msg - - l1 = pta.get_lnlikelihood(ps) - l2 = ptac.get_lnlikelihood(psc) - - # I don't know how to integrate l2 to match l1... - msg = 'Marginal and hierarchical likelihoods should be different.' - assert l1 != l2, msg + db = ptac.get_delay(params) + + msg = 'Implicit and explicit GP delays are different.' + assert np.allclose(da[0], db[0]), msg + assert np.allclose(da[1], db[1]), msg + + cpar = [p for p in ptac.params if 'coefficients' in p.name] + + def shouldfail(): + return cpar[0].get_logpdf(params) + + self.assertRaises(NotImplementedError, shouldfail) + + def test_fourier_red_noise(self): + """Test that implicit and explicit GP delays are the same.""" + # set up signal parameter + pl = utils.powerlaw(log10_A=parameter.Uniform(-18,-12), + gamma=parameter.Uniform(1,7)) + rn = gp_signals.FourierBasisGP(spectrum=pl, components=20) + rnm = rn(self.psr) + + rnc = gp_signals.FourierBasisGP(spectrum=pl, components=20, + coefficients=True) + rnmc = rnc(self.psr) + + # parameters + log10_A, gamma = -14.5, 4.33 + params = {'B1855+09_red_noise_log10_A': log10_A, + 'B1855+09_red_noise_gamma': gamma} + + # get the GP delays in two different ways + cf = np.random.randn(40) + d1 = np.dot(rnm.get_basis(params),cf) + + params.update({'B1855+09_red_noise_coefficients': cf}) + d2 = rnmc.get_delay(params) + + msg = 'Implicit and explicit GP delays are different.' + assert np.allclose(d1, d2), msg + + # np.array cast is needed because we get a KernelArray + phimat = np.array(rnm.get_phi(params)) + pr1 = (-0.5 * np.sum(cf*cf/phimat) - + 0.5 * np.sum(np.log(phimat)) - + 0.5 * len(phimat) * np.log(2*math.pi)) + + cpar = [p for p in rnmc.params if 'coefficients' in p.name][0] + pr2 = cpar.get_logpdf(params=params) + + msg = 'Implicit and explicit GP priors are different.' + assert np.allclose(pr1, pr2), msg + + def test_ecorr_backend(self): + """Test that ecorr-backend signal returns correct values.""" + # set up signal parameter + ecorr = parameter.Uniform(-10, -5) + selection = Selection(selections.by_backend) + ec = gp_signals.EcorrBasisModel(log10_ecorr=ecorr, + selection=selection) + ecm = ec(self.psr) + + ecc = gp_signals.EcorrBasisModel(log10_ecorr=ecorr, + selection=selection, + coefficients=True) + eccm = ecc(self.psr) + + # parameters + ecorrs = [-6.1, -6.2, -6.3, -6.4] + params = {'B1855+09_basis_ecorr_430_ASP_log10_ecorr': ecorrs[0], + 'B1855+09_basis_ecorr_430_PUPPI_log10_ecorr': ecorrs[1], + 'B1855+09_basis_ecorr_L-wide_ASP_log10_ecorr': ecorrs[2], + 'B1855+09_basis_ecorr_L-wide_PUPPI_log10_ecorr': + ecorrs[3]} + + fmat = ecm.get_basis(params) + cf = 1e-6 * np.random.randn(fmat.shape[1]) + d1 = np.dot(fmat, cf) + + for key in ecm._keys: + parname = 'B1855+09_basis_ecorr_' + key + '_coefficients' + params[parname] = cf[ecm._slices[key]] + d2 = eccm.get_delay(params) + + msg = 'Implicit and explicit ecorr-basis delays are different.' + assert np.allclose(d1, d2), msg + + def test_formalism(self): + # create marginalized model + ef = white_signals.MeasurementNoise( + efac=parameter.Uniform(0.1, 5.0)) + tm = gp_signals.TimingModel() + ec = gp_signals.EcorrBasisModel( + log10_ecorr=parameter.Uniform(-10, -5)) + pl = utils.powerlaw(log10_A=parameter.Uniform(-18,-12), + gamma=parameter.Uniform(1,7)) + rn = gp_signals.FourierBasisGP(spectrum=pl, components=10) + model = ef + tm + ec + rn + pta = signal_base.PTA([model(self.psr)]) + + # create hierarchical model + tmc = gp_signals.TimingModel(coefficients=True) + ecc = gp_signals.EcorrBasisModel( + log10_ecorr=parameter.Uniform(-10,-5), + coefficients=True) + rnc = gp_signals.FourierBasisGP(spectrum=pl, components=10, + coefficients=True) + modelc = ef + tmc + ecc + rnc + ptac = signal_base.PTA([modelc(self.psr)]) + + ps = {'B1855+09_efac': 1, + 'B1855+09_basis_ecorr_log10_ecorr': -6, + 'B1855+09_red_noise_log10_A': -14, + 'B1855+09_red_noise_gamma': 3} + psc = utils.get_coefficients(pta, ps) + + d1 = ptac.get_delay(psc)[0] + d2 = (np.dot(pta.pulsarmodels[0].signals[1].get_basis(ps), + psc['B1855+09_linear_timing_model_coefficients']) + + np.dot(pta.pulsarmodels[0].signals[2].get_basis(ps), + psc['B1855+09_basis_ecorr_coefficients']) + + np.dot(pta.pulsarmodels[0].signals[3].get_basis(ps), + psc['B1855+09_red_noise_coefficients'])) + + msg = 'Implicit and explicit PTA delays are different.' + assert np.allclose(d1, d2), msg + + l1 = pta.get_lnlikelihood(ps) + l2 = ptac.get_lnlikelihood(psc) + + # I don't know how to integrate l2 to match l1... + msg = 'Marginal and hierarchical likelihoods should be different.' + assert l1 != l2, msg class TestGPCoefficientsPint(TestGPCoefficients):