Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Simulate correlated DM noise in wideband TOAs #1868

Merged
merged 13 commits into from
Jan 17, 2025
Merged

Simulate correlated DM noise in wideband TOAs #1868

Show file tree
Hide file tree
Changes from 11 commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
99862c4
introduces_dm_errors
abhisrkckl Dec 5, 2024
90d6539
Simulate correlated DM noise for wideband TOAss
abhisrkckl Dec 5, 2024
45ff99f
fix for
abhisrkckl Dec 5, 2024
fdf0ba9
test_fake_toas
abhisrkckl Dec 5, 2024
fb94676
fix function call
abhisrkckl Dec 5, 2024
7bfb6f8
Merge branch 'master' into sim-dmgp
abhisrkckl Dec 6, 2024
0f37e15
cleanup
abhisrkckl Dec 6, 2024
29feb02
Merge branch 'nanograv:master' into sim-dmgp
abhisrkckl Dec 9, 2024
bffe712
Merge branch 'master' into sim-dmgp
abhisrkckl Dec 11, 2024
f6cbcd5
Merge branch 'nanograv:master' into sim-dmgp
abhisrkckl Dec 18, 2024
403433b
Merge branch 'master' into sim-dmgp
abhisrkckl Dec 19, 2024
ce48d47
Merge branch 'master' into sim-dmgp
abhisrkckl Jan 16, 2025
eaacbb4
test
abhisrkckl Jan 17, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGELOG-unreleased.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@ the released changes.
## Unreleased
### Changed
### Added
- Simulate correlated DM noise for wideband TOAs
### Fixed
### Removed

5 changes: 5 additions & 0 deletions src/pint/models/noise_model.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,9 @@


class NoiseComponent(Component):

introduces_dm_errors = False

def __init__(
self,
):
Expand Down Expand Up @@ -232,6 +235,7 @@ class ScaleDmError(NoiseComponent):
category = "scale_dm_error"

introduces_correlated_errors = False
introduces_dm_errors = True
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is the purpose of this?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Currently there is no way to distinguish between achromatic noise models and noise models that represent DM noise (such as PLDMNoise) that affects the DM part of the wideband likelihood. I added this attribute to distinguish such noise components. The timing model components that introduce a DM correction can already be distinguished by checking if they are derived from the Dispersion class.


def __init__(
self,
Expand Down Expand Up @@ -470,6 +474,7 @@ class PLDMNoise(NoiseComponent):
category = "pl_DM_noise"

introduces_correlated_errors = True
introduces_dm_errors = True
is_time_correlated = True

def __init__(
Expand Down
26 changes: 23 additions & 3 deletions src/pint/simulation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
from loguru import logger as log
from astropy import time

from pint.models.noise_model import NoiseComponent
from pint.types import time_like, file_like
import pint.residuals
import pint.toa
Expand Down Expand Up @@ -185,6 +186,7 @@ def update_fake_dms(
ts: pint.toa.TOAs,
dm_error: u.Quantity,
add_noise: bool,
add_correlated_noise: bool,
) -> pint.toa.TOAs:
"""Update simulated wideband DM information in TOAs.

Expand All @@ -209,6 +211,20 @@ def update_fake_dms(
if add_noise:
dms += scaled_dm_errors.to(pint.dmu) * np.random.randn(len(scaled_dm_errors))

if add_correlated_noise:
dm_noise = np.zeros(len(toas)) * pint.dmu
for noise_comp in model.NoiseComponent_list:
if (
noise_comp.introduces_correlated_errors
and noise_comp.introduces_dm_errors
):
U = noise_comp.get_noise_basis(toas)
b = noise_comp.get_noise_weights(toas)
delay = (U @ (b**0.5 * np.random.normal(size=len(b)))) << u.s
freqs = model.barycentric_radio_freq(toas)
dm_noise += (delay / pint.DMconst * freqs**2).to(pint.dmu)
dms += dm_noise

for f, dm in zip(toas.table["flags"], dms):
f["pp_dm"] = str(dm.to_value(pint.dmu))

Expand Down Expand Up @@ -338,7 +354,9 @@ def make_fake_toas_uniform(
)

if wideband:
ts = update_fake_dms(model, ts, wideband_dm_error, add_noise)
ts = update_fake_dms(
model, ts, wideband_dm_error, add_noise, add_correlated_noise
)

return make_fake_toas(
ts,
Expand Down Expand Up @@ -466,7 +484,9 @@ def make_fake_toas_fromMJDs(
)

if wideband:
ts = update_fake_dms(model, ts, wideband_dm_error, add_noise)
ts = update_fake_dms(
model, ts, wideband_dm_error, add_noise, add_correlated_noise
)

return make_fake_toas(
ts,
Expand Down Expand Up @@ -517,7 +537,7 @@ def make_fake_toas_fromtim(

if ts.is_wideband():
dm_errors = ts.get_dm_errors()
ts = update_fake_dms(model, ts, dm_errors, add_noise)
ts = update_fake_dms(model, ts, dm_errors, add_noise, add_correlated_noise)

return make_fake_toas(
ts,
Expand Down
54 changes: 54 additions & 0 deletions tests/test_fake_toas.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -511,3 +511,57 @@ def test_simulate_uniform_multifreq(multifreq):
multi_freqs_in_epoch=multifreq,
)
assert len(t) == ntoas


def test_simulate_wideband_dmgp():
par = """
PSRJ J0023+0923
RAJ 00:23:16.8790858 1 0.00002408141295805134
DECJ +09:23:23.86936 1 0.00082010713730773120
F0 327.84702062954047136 1 0.00000000000295205483
F1 -1.2278326306812866375e-15 1 3.8219244605614075223e-19
PEPOCH 56199.999797564144902
POSEPOCH 56199.999797564144902
DMEPOCH 56200
DM 14.327978186774068347 1 0.00006751663559857748
BINARY ELL1
PB 0.13879914244858396754 1 0.00000000003514075083
A1 0.034841158415224894973 1 0.00000012173038389247
TASC 56178.804891768506529 1 0.00000007765191894742
EPS1 1.6508830631753595232e-05 1 0.00000477568412215803
EPS2 3.9656838708709247373e-06 1 0.00000458951091435993
CLK TT(BIPM2015)
MODE 1
UNITS TDB
TIMEEPH FB90
DILATEFREQ N
PLANET_SHAPIRO N
CORRECT_TROPOSPHERE N
EPHEM DE436
TNRedAmp -13.3087
TNRedGam 1.5
TNRedC 14
TNDMAMP -12.2
TNDMGAM 3.5
TNDMC 15
"""

m = get_model(io.StringIO(par))
t = pint.simulation.make_fake_toas_uniform(
startMJD=54000,
endMJD=56000,
ntoas=200,
model=m,
add_noise=True,
wideband=True,
add_correlated_noise=True,
)

assert t.is_wideband()

# There is correlated noise in DMs
assert (
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So this is some sort of a chi^2 test? Is there no better way to check for correlated noise?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. A better test will be something like Ljung-Box test. Unfortunately, it is not available in scipy. We will need to add statsmodels as a dependency to use that test.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It looks like that test would not be very hard to implement, maybe include in PINT utils. Although it may not be needed just for this, it could be more broadly useful.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK... I have added the Ljung-Box test. For this I added statsmodels as a dependency for tests. It is not a dependency for the package itself since it is not used anywhere in the main code. #1776 may change that though.

sum(((t.get_dms() - m.total_dm(t)) / m.scaled_dm_uncertainty(t)) ** 2).si.value
/ len(t)
> 10
)
Loading