From cd03628aa7ec0ca7b89527ccda31be9480671a8d Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Tue, 26 Sep 2023 11:01:08 -0500 Subject: [PATCH 01/10] simrednoise --- src/pint/simulation.py | 39 ++++++++++++++++++++++++++++++++++----- 1 file changed, 34 insertions(+), 5 deletions(-) diff --git a/src/pint/simulation.py b/src/pint/simulation.py index 45b9eb5c8..c004a02bd 100644 --- a/src/pint/simulation.py +++ b/src/pint/simulation.py @@ -128,7 +128,7 @@ def get_fake_toa_clock_versions(model, include_bipm=False, include_gps=True): } -def make_fake_toas(ts, model, add_noise=False, name="fake"): +def make_fake_toas(ts, model, add_noise=False, add_correlated_noise=False, name="fake"): """Make toas from an array of times Can include alternating frequencies if fed an array of frequencies, @@ -156,6 +156,13 @@ def make_fake_toas(ts, model, add_noise=False, name="fake"): """ tsim = deepcopy(ts) zero_residuals(tsim, model) + + if add_correlated_noise: + U = model.noise_model_designmatrix(tsim) + b = model.noise_model_basis_weight(tsim) + corrn = (U @ (b**0.5 * np.random.randn(len(b)))) << u.s + tsim.adjust_TOAs(time.TimeDelta(corrn)) + if add_noise: # this function will include EFAC and EQUAD err = model.scaled_toa_uncertainty(tsim) * np.random.normal(size=len(tsim)) @@ -198,6 +205,7 @@ def make_fake_toas_uniform( obs="GBT", error=1 * u.us, add_noise=False, + add_correlated_noise=False, wideband=False, wideband_dm_error=1e-4 * pint.dmu, name="fake", @@ -301,7 +309,13 @@ def make_fake_toas_uniform( if wideband: ts = update_fake_dms(model, ts, wideband_dm_error, add_noise) - return make_fake_toas(ts, model=model, add_noise=add_noise, name=name) + return make_fake_toas( + ts, + model=model, + add_noise=add_noise, + add_correlated_noise=add_correlated_noise, + name=name, + ) def make_fake_toas_fromMJDs( @@ -311,6 +325,7 @@ def make_fake_toas_fromMJDs( obs="GBT", error=1 * u.us, add_noise=False, + add_correlated_noise=False, wideband=False, wideband_dm_error=1e-4 * pint.dmu, name="fake", @@ -396,10 +411,18 @@ def make_fake_toas_fromMJDs( if wideband: ts = update_fake_dms(model, ts, wideband_dm_error, add_noise) - return make_fake_toas(ts, model=model, add_noise=add_noise, name=name) + return make_fake_toas( + ts, + model=model, + add_noise=add_noise, + add_correlated_noise=add_correlated_noise, + name=name, + ) -def make_fake_toas_fromtim(timfile, model, add_noise=False, name="fake"): +def make_fake_toas_fromtim( + timfile, model, add_noise=False, add_correlated_noise=False, name="fake" +): """Make fake toas with the same times as an input tim file Can include alternating frequencies if fed an array of frequencies, @@ -432,7 +455,13 @@ def make_fake_toas_fromtim(timfile, model, add_noise=False, name="fake"): dm_errors = input_ts.get_dm_errors() ts = update_fake_dms(model, ts, dm_errors, add_noise) - return make_fake_toas(input_ts, model=model, add_noise=add_noise, name=name) + return make_fake_toas( + input_ts, + model=model, + add_noise=add_noise, + add_correlated_noise=add_correlated_noise, + name=name, + ) def calculate_random_models( From c381e568a32541e2726e746004d046b2ba9ae17f Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Tue, 26 Sep 2023 11:04:27 -0500 Subject: [PATCH 02/10] zim --- src/pint/scripts/zima.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/pint/scripts/zima.py b/src/pint/scripts/zima.py index 2f5f370a5..89641ee07 100755 --- a/src/pint/scripts/zima.py +++ b/src/pint/scripts/zima.py @@ -65,6 +65,12 @@ def main(argv=None): default=False, help="Actually add in random noise, or just populate the column", ) + parser.add_argument( + "--addcorrnoise", + action="store_true", + default=False, + help="Add in a correlated noise realization if it's present in the model", + ) parser.add_argument( "--wideband", action="store_true", @@ -127,13 +133,17 @@ def main(argv=None): freq=np.atleast_1d(args.freq) * u.MHz, fuzz=args.fuzzdays * u.d, add_noise=args.addnoise, + add_correlated_noise=args.addcorrnoise, wideband=args.wideband, wideband_dm_error=args.dmerror * pint.dmu, ) else: log.info(f"Reading initial TOAs from {args.inputtim}") ts = pint.simulation.make_fake_toas_fromtim( - args.inputtim, model=m, add_noise=args.addnoise + args.inputtim, + model=m, + add_noise=args.addnoise, + add_correlated_noise=args.addcorrnoise, ) # Write TOAs to a file From c6f8839f2d0d27aa6060f202f76f1c9d828fff45 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Tue, 26 Sep 2023 11:05:41 -0500 Subject: [PATCH 03/10] CHANGELOG --- CHANGELOG-unreleased.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 1bebb820a..2b67108aa 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -18,6 +18,7 @@ the released changes. - `pint.output.publish` module and `pintpublish` script for generating publication (LaTeX) output. - Added radial velocity methods for binary models - Added `DMWaveX` model (Fourier representation of DM noise) +- Simulate correlated noise using `pint.simulate` functions and the `zima` script. ### Fixed - Wave model `validate()` can correctly use PEPOCH to assign WAVEEPOCH parameter - Fixed RTD by specifying theme explicitly. From d5270f929bc9051f6e514bb9b04b065d958ea2cd Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Tue, 26 Sep 2023 13:47:09 -0500 Subject: [PATCH 04/10] -- --- tests/.gitignore | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 tests/.gitignore diff --git a/tests/.gitignore b/tests/.gitignore new file mode 100644 index 000000000..f1cf55bf1 --- /dev/null +++ b/tests/.gitignore @@ -0,0 +1,2 @@ +*.ipynb +chains/* From 07ebfbed66c9fbfedb57a4771a10ab0263d21c9e Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Tue, 26 Sep 2023 14:01:03 -0500 Subject: [PATCH 05/10] test_zima --- src/pint/data/examples/B1855+09_NANOGrav_9yv1.gls.par | 2 +- tests/datafile/B1855+09_NANOGrav_9yv1.gls.par | 2 +- tests/test_zima.py | 8 ++++++++ 3 files changed, 10 insertions(+), 2 deletions(-) diff --git a/src/pint/data/examples/B1855+09_NANOGrav_9yv1.gls.par b/src/pint/data/examples/B1855+09_NANOGrav_9yv1.gls.par index 2138184bd..645cb802a 100644 --- a/src/pint/data/examples/B1855+09_NANOGrav_9yv1.gls.par +++ b/src/pint/data/examples/B1855+09_NANOGrav_9yv1.gls.par @@ -453,7 +453,7 @@ FD3 1.07526915D-04 1 2.50177766D-05 SOLARN0 0.00 EPHEM DE421 ECL IERS2003 -CLK TT(BIPM) +CLK TT(BIPM2019) UNITS TDB TIMEEPH FB90 #T2CMETHOD TEMPO diff --git a/tests/datafile/B1855+09_NANOGrav_9yv1.gls.par b/tests/datafile/B1855+09_NANOGrav_9yv1.gls.par index 2138184bd..645cb802a 100644 --- a/tests/datafile/B1855+09_NANOGrav_9yv1.gls.par +++ b/tests/datafile/B1855+09_NANOGrav_9yv1.gls.par @@ -453,7 +453,7 @@ FD3 1.07526915D-04 1 2.50177766D-05 SOLARN0 0.00 EPHEM DE421 ECL IERS2003 -CLK TT(BIPM) +CLK TT(BIPM2019) UNITS TDB TIMEEPH FB90 #T2CMETHOD TEMPO diff --git a/tests/test_zima.py b/tests/test_zima.py index d05664669..dc6a5f387 100644 --- a/tests/test_zima.py +++ b/tests/test_zima.py @@ -92,3 +92,11 @@ def test_zima_fuzzdays(tmp_path): lines = sys.stdout.getvalue() finally: sys.stdout = saved_stdout + + +def test_zima_corrnoise(tmp_path): + parfile = datadir / "B1855+09_NANOGrav_9yv1.gls.par" + output_timfile = tmp_path / "fake_testzima.tim" + cmd = f"--addnoise --addcorrnoise {parfile} {output_timfile}" + zima.main(cmd.split()) + assert os.path.isfile(output_timfile) From 24d8bafc58c5fab193f91674bb772296a951b15d Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Tue, 26 Sep 2023 14:06:26 -0500 Subject: [PATCH 06/10] docs --- src/pint/simulation.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/pint/simulation.py b/src/pint/simulation.py index c004a02bd..e57bd50b0 100644 --- a/src/pint/simulation.py +++ b/src/pint/simulation.py @@ -142,6 +142,8 @@ def make_fake_toas(ts, model, add_noise=False, add_correlated_noise=False, name= current model add_noise : bool, optional Add noise to the TOAs (otherwise `error` just populates the column) + add_correlated_noise : bool, optional + Add correlated noise to the TOAs if it's present in the timing mode. name : str, optional Name for the TOAs (goes into the flags) @@ -237,6 +239,8 @@ def make_fake_toas_uniform( uncertainty to attach to each TOA add_noise : bool, optional Add noise to the TOAs (otherwise `error` just populates the column) + add_correlated_noise : bool, optional + Add correlated noise to the TOAs if it's present in the timing mode. wideband : bool, optional Whether to include wideband DM information with each TOA; default is not to include any wideband DM information. If True, the DM associated @@ -351,6 +355,8 @@ def make_fake_toas_fromMJDs( uncertainty to attach to each TOA add_noise : bool, optional Add noise to the TOAs (otherwise `error` just populates the column) + add_correlated_noise : bool, optional + Add correlated noise to the TOAs if it's present in the timing mode. wideband : astropy.units.Quantity, optional Whether to include wideband DM values with each TOA; default is not to include any DM information @@ -436,6 +442,8 @@ def make_fake_toas_fromtim( current model add_noise : bool, optional Add noise to the TOAs (otherwise `error` just populates the column) + add_correlated_noise : bool, optional + Add correlated noise to the TOAs if it's present in the timing mode. name : str, optional Name for the TOAs (goes into the flags) From 62f7557399346cf0db664a12d69d7700182f47d2 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 28 Sep 2023 16:09:31 -0500 Subject: [PATCH 07/10] remove gitignore --- src/pint/simulation.py | 2 +- tests/.gitignore | 2 -- 2 files changed, 1 insertion(+), 3 deletions(-) delete mode 100644 tests/.gitignore diff --git a/src/pint/simulation.py b/src/pint/simulation.py index e57bd50b0..92a9c71d1 100644 --- a/src/pint/simulation.py +++ b/src/pint/simulation.py @@ -162,7 +162,7 @@ def make_fake_toas(ts, model, add_noise=False, add_correlated_noise=False, name= if add_correlated_noise: U = model.noise_model_designmatrix(tsim) b = model.noise_model_basis_weight(tsim) - corrn = (U @ (b**0.5 * np.random.randn(len(b)))) << u.s + corrn = (U @ (b**0.5 * np.random.normal(size=len(b)))) << u.s tsim.adjust_TOAs(time.TimeDelta(corrn)) if add_noise: diff --git a/tests/.gitignore b/tests/.gitignore deleted file mode 100644 index f1cf55bf1..000000000 --- a/tests/.gitignore +++ /dev/null @@ -1,2 +0,0 @@ -*.ipynb -chains/* From 0bb171d264f17a1ae64f7c1aa7c3691d6272b57f Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 28 Sep 2023 17:17:02 -0500 Subject: [PATCH 08/10] test_simulate_corrnoise --- tests/test_zima.py | 44 ++++++++++++++++++++++++++++++++++++++------ 1 file changed, 38 insertions(+), 6 deletions(-) diff --git a/tests/test_zima.py b/tests/test_zima.py index dc6a5f387..bafbb3206 100644 --- a/tests/test_zima.py +++ b/tests/test_zima.py @@ -8,8 +8,10 @@ import matplotlib import pint.scripts.zima as zima -from pint.models import get_model_and_toas +from pint.models import get_model_and_toas, get_model +from pint.simulation import make_fake_toas_uniform from pint.residuals import Residuals +from pint.fitter import DownhillGLSFitter @pytest.mark.parametrize("addnoise", ["", "--addnoise"]) @@ -94,9 +96,39 @@ def test_zima_fuzzdays(tmp_path): sys.stdout = saved_stdout -def test_zima_corrnoise(tmp_path): +def test_simulate_corrnoise(tmp_path): parfile = datadir / "B1855+09_NANOGrav_9yv1.gls.par" - output_timfile = tmp_path / "fake_testzima.tim" - cmd = f"--addnoise --addcorrnoise {parfile} {output_timfile}" - zima.main(cmd.split()) - assert os.path.isfile(output_timfile) + + m = get_model(parfile) + + # Simulated TOAs won't have the correct flags for some of these to work. + m.remove_component("ScaleToaError") + m.remove_component("EcorrNoise") + m.remove_component("DispersionDMX") + m.remove_component("PhaseJump") + m.remove_component("FD") + m.PLANET_SHAPIRO.value = False + + t = make_fake_toas_uniform( + m.START.value, + m.FINISH.value, + 1000, + m, + add_noise=True, + add_correlated_noise=True, + ) + + # Check if the created TOAs can be whitened using + # the original timing model. This won't work if the + # noise is not realized correctly. + ftr = DownhillGLSFitter(t, m) + ftr.fit_toas() + rc = sum(ftr.resids.noise_resids.values()) + r = ftr.resids.time_resids + rw = r - rc + sigma = ftr.resids.get_data_error() + + # This should be independent and standard-normal distributed. + x = (rw / sigma).to_value("") + assert np.isclose(np.std(x), 1, atol=0.2) + assert np.isclose(np.mean(x), 0, atol=0.01) From 2e17704fdc18937b38c2dc33a55120285b7cd776 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Mon, 9 Oct 2023 16:44:56 -0500 Subject: [PATCH 09/10] fix wavex_setup and dmwavex_setup --- src/pint/utils.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/pint/utils.py b/src/pint/utils.py index c4410618a..9fcdb31dc 100644 --- a/src/pint/utils.py +++ b/src/pint/utils.py @@ -1328,8 +1328,9 @@ def wavex_setup(model, T_span, freqs=None, n_freqs=None): "Both freqs and n_freqs are specified. Only one or the other should be used" ) - if n_freqs <= 0: + if n_freqs is not None and n_freqs <= 0: raise ValueError("Must use a non-zero number of wave frequencies") + model.add_component(WaveX()) if isinstance(T_span, u.quantity.Quantity): T_span.to(u.d) @@ -1356,11 +1357,11 @@ def wavex_setup(model, T_span, freqs=None, n_freqs=None): if n_freqs is not None: if n_freqs == 1: - wave_freq = 2.0 * np.pi / T_span + wave_freq = 1 / T_span model.WXFREQ_0001.quantity = wave_freq else: wave_numbers = np.arange(1, n_freqs + 1) - wave_freqs = 2.0 * np.pi * wave_numbers / T_span + wave_freqs = wave_numbers / T_span model.WXFREQ_0001.quantity = wave_freqs[0] model.components["WaveX"].add_wavex_components(wave_freqs[1:]) return model.components["WaveX"].get_indices() @@ -1408,8 +1409,9 @@ def dmwavex_setup(model, T_span, freqs=None, n_freqs=None): "Both freqs and n_freqs are specified. Only one or the other should be used" ) - if n_freqs <= 0: + if n_freqs is not None and n_freqs <= 0: raise ValueError("Must use a non-zero number of wave frequencies") + model.add_component(DMWaveX()) if isinstance(T_span, u.quantity.Quantity): T_span.to(u.d) @@ -1436,11 +1438,11 @@ def dmwavex_setup(model, T_span, freqs=None, n_freqs=None): if n_freqs is not None: if n_freqs == 1: - wave_freq = 2.0 * np.pi / T_span + wave_freq = 1 / T_span model.DMWXFREQ_0001.quantity = wave_freq else: wave_numbers = np.arange(1, n_freqs + 1) - wave_freqs = 2.0 * np.pi * wave_numbers / T_span + wave_freqs = wave_numbers / T_span model.DMWXFREQ_0001.quantity = wave_freqs[0] model.components["DMWaveX"].add_dmwavex_components(wave_freqs[1:]) return model.components["DMWaveX"].get_indices() From d7f93e426606974c1c1f19b02bea8350b23b6809 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Tue, 10 Oct 2023 10:03:14 -0500 Subject: [PATCH 10/10] Update AUTHORS.rst --- AUTHORS.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/AUTHORS.rst b/AUTHORS.rst index 07bbac26c..de5bf4441 100644 --- a/AUTHORS.rst +++ b/AUTHORS.rst @@ -41,7 +41,7 @@ Active developers are indicated by (*). Authors of the PINT paper are indicated * Sasha Levina * Nikhil Mahajan * Alex McEwen -* Patrick O'Neill +* Patrick O'Neill (*) * Tim Pennucci * Camryn Phillips (#) * Matt Pitkin