From d6611a3832022da60431bbd9859f845eeb663dfd Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Thu, 31 Aug 2023 09:17:20 -0500 Subject: [PATCH 1/8] include CHI2 in par output --- CHANGELOG-unreleased.md | 1 + src/pint/fitter.py | 1 + tests/test_chi2inpar.py | 47 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 49 insertions(+) create mode 100644 tests/test_chi2inpar.py diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index ac9bc2495..b3cafbba6 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -12,6 +12,7 @@ the released changes. - Moved design matrix normalization code from `pint.fitter` to the new `pint.utils.normalize_designmatrix()` function. - Made `Residuals` independent of `GLSFitter` (GLS chi2 is now computed using the new function `Residuals._calc_gls_chi2()`). ### Added +- CHI2 now in par files ### Fixed - Fixed RTD by specifying theme explicitly. - `.value()` now works for pairParameters diff --git a/src/pint/fitter.py b/src/pint/fitter.py index 61bcd9532..bd754d467 100644 --- a/src/pint/fitter.py +++ b/src/pint/fitter.py @@ -698,6 +698,7 @@ def update_model(self, chi2=None): if self.toas.clock_corr_info["include_bipm"] else "TT(TAI)" ) + self.model.CHI2.value = self.resids.chi2 def reset_model(self): """Reset the current model to the initial model.""" diff --git a/tests/test_chi2inpar.py b/tests/test_chi2inpar.py new file mode 100644 index 000000000..ec8e908fe --- /dev/null +++ b/tests/test_chi2inpar.py @@ -0,0 +1,47 @@ +import pytest +from pint.models import get_model_and_toas, get_model +import pint.fitter +from pinttestdata import datadir +import os +import io + + +@pytest.mark.parametrize( + "ft", + [ + pint.fitter.WLSFitter, + pint.fitter.GLSFitter, + pint.fitter.DownhillWLSFitter, + pint.fitter.DownhillGLSFitter, + ], +) +def test_nbfitters(ft): + m, t = get_model_and_toas( + os.path.join(datadir, "NGC6440E.par"), os.path.join(datadir, "NGC6440E.tim") + ) + f = ft(t, m) + f.fit_toas() + matched_line = [ + line for line in f.model.as_parfile().split("\n") if line.startswith("CHI2") + ] + assert matched_line + assert float(matched_line[0].split()[-1]) > 0 + m2 = get_model(io.StringIO(f.model.as_parfile())) + + +@pytest.mark.parametrize( + "ft", [pint.fitter.WidebandTOAFitter, pint.fitter.WidebandDownhillFitter] +) +def test_wbfitters(ft): + m, t = get_model_and_toas( + "tests/datafile/B1855+09_NANOGrav_12yv3.wb.gls.par", + "tests/datafile/B1855+09_NANOGrav_12yv3.wb.tim", + ) + f = ft(t, m) + f.fit_toas() + matched_line = [ + line for line in f.model.as_parfile().split("\n") if line.startswith("CHI2") + ] + assert matched_line + assert float(matched_line[0].split()[-1]) > 0 + m2 = get_model(io.StringIO(f.model.as_parfile())) From a63ffa08d6e5071f3ecb05fc5a35b470f8b0fec1 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Thu, 31 Aug 2023 09:39:55 -0500 Subject: [PATCH 2/8] made tests faster --- tests/test_chi2inpar.py | 35 ++++++++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/tests/test_chi2inpar.py b/tests/test_chi2inpar.py index ec8e908fe..fa389d6d7 100644 --- a/tests/test_chi2inpar.py +++ b/tests/test_chi2inpar.py @@ -1,11 +1,33 @@ +from astropy import units as u import pytest from pint.models import get_model_and_toas, get_model import pint.fitter +import pint.simulation from pinttestdata import datadir import os import io +@pytest.fixture +def wb(): + m = get_model(os.path.join(datadir, "NGC6440E.par")) + t = pint.simulation.make_fake_toas_uniform( + 55000, 58000, 20, model=m, freq=1400 * u.MHz, wideband=True + ) + + return m, t + + +@pytest.fixture +def nb(): + m = get_model(os.path.join(datadir, "NGC6440E.par")) + t = pint.simulation.make_fake_toas_uniform( + 55000, 58000, 20, model=m, freq=1400 * u.MHz, wideband=False + ) + + return m, t + + @pytest.mark.parametrize( "ft", [ @@ -15,10 +37,8 @@ pint.fitter.DownhillGLSFitter, ], ) -def test_nbfitters(ft): - m, t = get_model_and_toas( - os.path.join(datadir, "NGC6440E.par"), os.path.join(datadir, "NGC6440E.tim") - ) +def test_nbfitters(ft, nb): + m, t = nb f = ft(t, m) f.fit_toas() matched_line = [ @@ -32,11 +52,8 @@ def test_nbfitters(ft): @pytest.mark.parametrize( "ft", [pint.fitter.WidebandTOAFitter, pint.fitter.WidebandDownhillFitter] ) -def test_wbfitters(ft): - m, t = get_model_and_toas( - "tests/datafile/B1855+09_NANOGrav_12yv3.wb.gls.par", - "tests/datafile/B1855+09_NANOGrav_12yv3.wb.tim", - ) +def test_wbfitters(ft, wb): + m, t = wb f = ft(t, m) f.fit_toas() matched_line = [ From 99faa714fe24522c2c8ca2aaf92852d969424324 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Thu, 31 Aug 2023 09:58:57 -0500 Subject: [PATCH 3/8] removed infinite loop --- src/pint/fitter.py | 3 ++- tests/test_chi2inpar.py | 10 ++++++++-- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/src/pint/fitter.py b/src/pint/fitter.py index bd754d467..f6aa9755d 100644 --- a/src/pint/fitter.py +++ b/src/pint/fitter.py @@ -698,7 +698,8 @@ def update_model(self, chi2=None): if self.toas.clock_corr_info["include_bipm"] else "TT(TAI)" ) - self.model.CHI2.value = self.resids.chi2 + if chi2 is not None: + self.model.CHI2.value = chi2 def reset_model(self): """Reset the current model to the initial model.""" diff --git a/tests/test_chi2inpar.py b/tests/test_chi2inpar.py index fa389d6d7..69263ebfd 100644 --- a/tests/test_chi2inpar.py +++ b/tests/test_chi2inpar.py @@ -12,7 +12,7 @@ def wb(): m = get_model(os.path.join(datadir, "NGC6440E.par")) t = pint.simulation.make_fake_toas_uniform( - 55000, 58000, 20, model=m, freq=1400 * u.MHz, wideband=True + 55000, 58000, 20, model=m, freq=1400 * u.MHz, wideband=True, add_noise=True ) return m, t @@ -22,7 +22,13 @@ def wb(): def nb(): m = get_model(os.path.join(datadir, "NGC6440E.par")) t = pint.simulation.make_fake_toas_uniform( - 55000, 58000, 20, model=m, freq=1400 * u.MHz, wideband=False + 55000, + 58000, + 20, + model=m, + freq=1400 * u.MHz, + wideband=False, + add_noise=True, ) return m, t From 83d61a1fe50576f69bbcdec0ba3b759356612144 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Thu, 31 Aug 2023 14:29:51 -0500 Subject: [PATCH 4/8] added CHI2R, DMRES, TRES --- CHANGELOG-unreleased.md | 2 +- src/pint/fitter.py | 7 +++++++ src/pint/models/timing_model.py | 31 +++++++++++++++++++++++++++---- tests/test_chi2inpar.py | 23 ++++++++++++++--------- 4 files changed, 49 insertions(+), 14 deletions(-) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index b3cafbba6..2bcb21210 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -12,7 +12,7 @@ the released changes. - Moved design matrix normalization code from `pint.fitter` to the new `pint.utils.normalize_designmatrix()` function. - Made `Residuals` independent of `GLSFitter` (GLS chi2 is now computed using the new function `Residuals._calc_gls_chi2()`). ### Added -- CHI2 now in par files +- CHI2, CHI2R, TRES, DMRES now in par files ### Fixed - Fixed RTD by specifying theme explicitly. - `.value()` now works for pairParameters diff --git a/src/pint/fitter.py b/src/pint/fitter.py index f6aa9755d..823170f14 100644 --- a/src/pint/fitter.py +++ b/src/pint/fitter.py @@ -699,7 +699,14 @@ def update_model(self, chi2=None): else "TT(TAI)" ) if chi2 is not None: + # assume a fit has been done self.model.CHI2.value = chi2 + self.model.CHI2R.value = chi2 / self.resids.dof + if not self.is_wideband: + self.model.TRES.quantity = self.resids.rms_weighted() + else: + self.model.TRES.quantity = self.resids.rms_weighted()["toa"] + self.model.DMRES.quantity = self.resids.rms_weighted()["dm"] def reset_model(self): """Reset the current model to the initial model.""" diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 57eba15fa..e07bb2430 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -89,13 +89,13 @@ # # Comparisons with keywords in par file lines is done in a case insensitive way. ignore_params = { - "TRES", + # "TRES", "TZRMJD", "TZRFRQ", "TZRSITE", "NITS", "IBOOT", - "CHI2R", + # "CHI2R", "MODE", "PLANET_SHAPIRO2", } @@ -338,13 +338,36 @@ def __init__(self, name="", components=[]): self.add_param_from_top( floatParameter( name="CHI2", - value=0.0, units="", description="Chi-squared value obtained during fitting", ), "", ) + self.add_param_from_top( + floatParameter( + name="CHI2R", + units="", + description="Reduced chi-squared value obtained during fitting", + ), + "", + ) + self.add_param_from_top( + floatParameter( + name="TRES", + units=u.us, + description="TOA residual after fitting", + ), + "", + ) + self.add_param_from_top( + floatParameter( + name="DMRES", + units=u.pc / u.cm**3, + description="DM residual after fitting (wideband only)", + ), + "", + ) for cp in components: self.add_component(cp, setup=False, validate=False) @@ -2203,7 +2226,7 @@ def compare( else: value2[pn] = str(otherpar.value) if otherpar.value != par.value: - if par.name in ["START", "FINISH", "CHI2", "NTOA"]: + if par.name in ["START", "FINISH", "CHI2", "CHI2R", "NTOA"]: if verbosity == "max": log.info( "Parameter %s has changed between these models" diff --git a/tests/test_chi2inpar.py b/tests/test_chi2inpar.py index 69263ebfd..f7862a7cd 100644 --- a/tests/test_chi2inpar.py +++ b/tests/test_chi2inpar.py @@ -47,11 +47,15 @@ def test_nbfitters(ft, nb): m, t = nb f = ft(t, m) f.fit_toas() - matched_line = [ - line for line in f.model.as_parfile().split("\n") if line.startswith("CHI2") + for p in ["CHI2", "CHI2R", "TRES"]: + matched_line = [ + line for line in f.model.as_parfile().split("\n") if line.startswith(p) + ] + assert matched_line + assert float(matched_line[0].split()[-1]) > 0 + assert not [ + line for line in f.model.as_parfile().split("\n") if line.startswith("DMRES") ] - assert matched_line - assert float(matched_line[0].split()[-1]) > 0 m2 = get_model(io.StringIO(f.model.as_parfile())) @@ -62,9 +66,10 @@ def test_wbfitters(ft, wb): m, t = wb f = ft(t, m) f.fit_toas() - matched_line = [ - line for line in f.model.as_parfile().split("\n") if line.startswith("CHI2") - ] - assert matched_line - assert float(matched_line[0].split()[-1]) > 0 + for p in ["CHI2", "CHI2R", "TRES", "DMRES"]: + matched_line = [ + line for line in f.model.as_parfile().split("\n") if line.startswith(p) + ] + assert matched_line + assert float(matched_line[0].split()[-1]) > 0 m2 = get_model(io.StringIO(f.model.as_parfile())) From 954ef832c6ab39a0a0c607f2d81e63d8b088657d Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Thu, 31 Aug 2023 18:10:14 -0500 Subject: [PATCH 5/8] fix summary print --- src/pint/fitter.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/src/pint/fitter.py b/src/pint/fitter.py index 823170f14..2e6a31a9b 100644 --- a/src/pint/fitter.py +++ b/src/pint/fitter.py @@ -415,17 +415,22 @@ def get_summary(self, nodmx=False): pn, prefitpar.str_quantity(prefitpar.value), "", par.units ) elif par.frozen: - if ( - par.name == "START" - and prefitpar.value is None - or par.name != "START" - and par.name == "FINISH" + if par.name in ["START", "FINISH"] and prefitpar.value is None: + s += ("{:" + spacingName + "s} {:20s} {:28g} {} \n").format( + pn, " ", par.value, par.units + ) + elif par.name in ["START", "FINISH"]: + s += ("{:" + spacingName + "s} {:20g} {:28g} {} \n").format( + pn, prefitpar.value, par.value, par.units + ) + elif ( + par.name in ["CHI2", "CHI2R", "TRES", "DMRES"] and prefitpar.value is None ): s += ("{:" + spacingName + "s} {:20s} {:28g} {} \n").format( pn, " ", par.value, par.units ) - elif par.name in ["START", "FINISH"]: + elif par.name in ["CHI2", "CHI2R", "TRES", "DMRES"]: s += ("{:" + spacingName + "s} {:20g} {:28g} {} \n").format( pn, prefitpar.value, par.value, par.units ) @@ -433,6 +438,7 @@ def get_summary(self, nodmx=False): s += ("{:" + spacingName + "s} {:20g} {:28s} {} \n").format( pn, prefitpar.value, "", par.units ) + else: # s += "{:14s} {:20g} {:20g} {:20.2g} {} \n".format( # pn, From c92de4d0ea39fa1e78dc5621910dce129a2b4001 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Thu, 31 Aug 2023 21:18:48 -0500 Subject: [PATCH 6/8] test requires fit; no CHI2R if not PINT format --- src/pint/models/parameter.py | 5 ++++- tests/test_parfile_writing_format.py | 1 + 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/pint/models/parameter.py b/src/pint/models/parameter.py index 016687563..78178e176 100644 --- a/src/pint/models/parameter.py +++ b/src/pint/models/parameter.py @@ -463,9 +463,12 @@ def as_parfile_line(self, format="pint"): name = self.name if self.use_alias is None else self.use_alias # special cases for parameter names that change depending on format - if self.name == "CHI2" and format.lower() != "pint": + if self.name in ["CHI2", "CHI2R"] and format.lower() != "pint": # no CHI2 for TEMPO/TEMPO2 return "" + if self.name in ["TRES", "DMRES"] and format.lower() not in ["pint", "tempo"]: + # TRES/DMRES only for PINT, TEMPO + return "" elif self.name == "SWM" and format.lower() != "pint": # no SWM for TEMPO/TEMPO2 return "" diff --git a/tests/test_parfile_writing_format.py b/tests/test_parfile_writing_format.py index f2795a1d5..4941c768c 100644 --- a/tests/test_parfile_writing_format.py +++ b/tests/test_parfile_writing_format.py @@ -29,6 +29,7 @@ def test_CHI2(): ) f = fitter.WLSFitter(toas=t, model=m) + f.fit_toas() assert "CHI2" in f.model.as_parfile() assert "CHI2" not in f.model.as_parfile(format="tempo2") assert "CHI2" not in f.model.as_parfile(format="tempo") From 7146f36a056f6876dd6e81a0e96b4c526e0668ef Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Thu, 14 Sep 2023 12:28:29 -0500 Subject: [PATCH 7/8] CHI2/CHI2R/TRES in all output parfile formats --- src/pint/models/parameter.py | 7 ++----- tests/test_parfile_writing_format.py | 4 ++-- 2 files changed, 4 insertions(+), 7 deletions(-) diff --git a/src/pint/models/parameter.py b/src/pint/models/parameter.py index 78178e176..eb66f5f61 100644 --- a/src/pint/models/parameter.py +++ b/src/pint/models/parameter.py @@ -463,11 +463,8 @@ def as_parfile_line(self, format="pint"): name = self.name if self.use_alias is None else self.use_alias # special cases for parameter names that change depending on format - if self.name in ["CHI2", "CHI2R"] and format.lower() != "pint": - # no CHI2 for TEMPO/TEMPO2 - return "" - if self.name in ["TRES", "DMRES"] and format.lower() not in ["pint", "tempo"]: - # TRES/DMRES only for PINT, TEMPO + if self.name in ["DMRES"] and format.lower() not in ["pint"]: + # DMRES only for PINT return "" elif self.name == "SWM" and format.lower() != "pint": # no SWM for TEMPO/TEMPO2 diff --git a/tests/test_parfile_writing_format.py b/tests/test_parfile_writing_format.py index 4941c768c..87b95d63a 100644 --- a/tests/test_parfile_writing_format.py +++ b/tests/test_parfile_writing_format.py @@ -31,8 +31,8 @@ def test_CHI2(): f = fitter.WLSFitter(toas=t, model=m) f.fit_toas() assert "CHI2" in f.model.as_parfile() - assert "CHI2" not in f.model.as_parfile(format="tempo2") - assert "CHI2" not in f.model.as_parfile(format="tempo") + assert "CHI2" in f.model.as_parfile(format="tempo2") + assert "CHI2" in f.model.as_parfile(format="tempo") def test_T2CMETHOD(): From 813acfae02afd6aa2a076e197b2d08706c6ae19b Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Mon, 25 Sep 2023 09:11:52 -0500 Subject: [PATCH 8/8] updated changelog --- CHANGELOG-unreleased.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index b68c184f3..3b754e927 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -13,7 +13,7 @@ the released changes. - Moved design matrix normalization code from `pint.fitter` to the new `pint.utils.normalize_designmatrix()` function. - Made `Residuals` independent of `GLSFitter` (GLS chi2 is now computed using the new function `Residuals._calc_gls_chi2()`). ### Added -- CHI2, CHI2R, TRES, DMRES now in par files +- CHI2, CHI2R, TRES, DMRES now in postfit par files - Added `WaveX` model as a `DelayComponent` with Fourier amplitudes as fitted parameters - `Parameter.as_latex` method for latex representation of a parameter. - `pint.output.publish` module and `pintpublish` script for generating publication (LaTeX) output.