Merge pull request #900 from aarchiba/fitter_update_model

Fitter now sets informational parameters

CHANGELOG.md

@@ -39,6 +39,7 @@ and this project, at least loosely, adheres to [Semantic Versioning](https://sem
 - Added support for wideband-TOA fitting (Pennucci 2019).
 - Added START and FINISH parameters as MJDParameters to timing_model. They are now modified after a fit and are displayed with a model's .par file output.
 - Added solar_angle calculation (PR #892)
+- Added parameters to TimingModel to support TEMPO/TEMPO2 compatible par files. (PR #900)
 - Added position vectors to Neptune (PR #901)
 - Added checking for TOAs in DMX bins and other similar parameters, if free (PR #874)
 - Added PiecewiseSpindown model, for spindown correction lasting for a given MJD range
@@ -51,6 +52,7 @@ and this project, at least loosely, adheres to [Semantic Versioning](https://sem
 - Changed requirements to astropy>=4.0
 - get_model can now read from file-like, including StringIO, objects (handy for testing) (PR #871)
 - WidebandDMResiduals now support access to their parts through .toa and .dm attributes (PR #861)
+- Fitters now update the fitted model to record things like EPHEM used for the TOAs (PR #900)
 - EFACs and EQUADs can be set independently from each other now. (PR #890)
 - WLSFitter and GLSFitter and WidebandTOAFitter can now report degenerate parameter combinations (PR #874)
 - Raise an exception if the DDK model is provided with SINI (PR #864)

src/pint/fitter.py

@@ -207,7 +207,10 @@ def get_summary(self, nodmx=False):
             if par.value is not None:
                 if isinstance(par, strParameter):
                     s += ("{:" + spacingName + "s} {:>20s} {:28s} {}\n").format(
-                        pn, prefitpar.value, "", par.units
+                        pn,
+                        prefitpar.value if prefitpar.value is not None else "",
+                        par.value,
+                        par.units,
                     )
                 elif isinstance(par, AngleParameter):
                     # Add special handling here to put uncertainty into arcsec
@@ -458,6 +461,25 @@ def plot(self):
         ax.grid(True)
         plt.show()
 
+    def update_model(self, chi2=None):
+        """Update the model to reflect fit results and TOA properties.
+
+        This is called by ``fit_toas`` to ensure that parameters like
+        ``START``, ``FINISH``, ``EPHEM``, and ``DMDATA`` are set in the model
+        to reflect the TOAs in actual use.
+        """
+        self.model.START.value = self.toas.first_MJD
+        self.model.FINISH.value = self.toas.last_MJD
+        self.model.NTOA.value = len(self.toas)
+        self.model.EPHEM.value = self.toas.ephem
+        self.model.DMDATA.value = hasattr(self.resids, "dm")
+        if not self.toas.clock_corr_info["include_bipm"]:
+            self.model.CLOCK.value = "TT(TAI)"
+        else:
+            self.model.CLOCK.value = f"TT({self.toas.clock_corr_info['bipm_version']})"
+        if chi2 is not None:
+            self.model.CHI2.value = chi2
+
     def reset_model(self):
         """Reset the current model to the initial model."""
         self.model = copy.deepcopy(self.model_init)
@@ -848,9 +870,7 @@ def fit_toas(self, maxiter=20):
         # necessarily the one that yields the best fit
         self.minimize_func(np.atleast_1d(self.fitresult.x), *list(fitp.keys()))
 
-        # Update START/FINISH params
-        self.model.START.value = self.toas.first_MJD
-        self.model.FINISH.value = self.toas.last_MJD
+        self.update_model(self.resids.chi2)
 
         return self.resids.chi2
 
@@ -989,9 +1009,7 @@ def fit_toas(self, maxiter=1, threshold=None):
             # Update Uncertainties
             self.set_param_uncertainties(fitperrs)
 
-        # Update START/FINISH params
-        self.model.START.value = self.toas.first_MJD
-        self.model.FINISH.value = self.toas.last_MJD
+        self.update_model(chi2)
 
         return chi2
 
@@ -1177,9 +1195,7 @@ def fit_toas(self, maxiter=1, threshold=0, full_cov=False):
                     noise_resids[comp] = np.dot(M[:, p0:p1], xhat[p0:p1]) * u.s
                 self.resids.noise_resids = noise_resids
 
-        # Update START/FINISH params
-        self.model.START.value = self.toas.first_MJD
-        self.model.FINISH.value = self.toas.last_MJD
+        self.update_model(chi2)
 
         return chi2
 
@@ -1510,8 +1526,6 @@ def fit_toas(self, maxiter=1, threshold=0, full_cov=False):
                     noise_resids[comp] = np.dot(M[:, p0:p1], xhat[p0:p1]) * u.s
                 self.resids.noise_resids = noise_resids
 
-        # Update START/FINISH params
-        self.model.START.value = self.toas.first_MJD
-        self.model.FINISH.value = self.toas.last_MJD
+        self.update_model(chi2)
 
         return chi2

src/pint/models/timing_model.py

@@ -9,6 +9,7 @@
 import inspect
 from collections import OrderedDict, defaultdict
 from functools import wraps
+from warnings import warn
 
 import astropy.time as time
 import astropy.units as u
@@ -21,6 +22,8 @@
 from pint.models.parameter import (
     AngleParameter,
     MJDParameter,
+    boolParameter,
+    floatParameter,
     Parameter,
     maskParameter,
     strParameter,
@@ -40,14 +43,6 @@
 # Comparisons with keywords in par file lines is done in a case insensitive way.
 ignore_params = set(
     [
-        "START",
-        "FINISH",
-        "EPHVER",
-        "UNITS",
-        "TIMEEPH",
-        "T2CMETHOD",
-        "DILATEFREQ",
-        "NTOA",
         "TRES",
         "TZRMJD",
         "TZRFRQ",
@@ -57,7 +52,6 @@
         "BINARY",
         "CHI2R",
         "MODE",
-        "INFO",
         "PLANET_SHAPIRO2",
         #    'NE_SW', 'NE_SW2',
     ]
@@ -228,6 +222,62 @@ def __init__(self, name="", components=[]):
         self.add_param_from_top(
             MJDParameter(name="FINISH", description="End MJD for fitting"), ""
         )
+        self.add_param_from_top(
+            strParameter(
+                name="INFO",
+                description="Tells TEMPO to write some extra information about frontend/backend combinations; -f is recommended",
+            ),
+            "",
+        )
+        self.add_param_from_top(
+            strParameter(
+                name="TIMEEPH",
+                description="Time ephemeris to use for TDB conversion; for PINT, always FB90",
+            ),
+            "",
+        )
+        self.add_param_from_top(
+            strParameter(
+                name="T2CMETHOD",
+                description="Method for transforming from terrestrial to celestial frame (IAU2000B/TEMPO; PINT only supports ????)",
+            ),
+            "",
+        )
+        self.add_param_from_top(
+            boolParameter(
+                name="DILATEFREQ",
+                value=False,
+                description="Whether or not TEMPO2 should apply gravitational redshift and time dilation to obseerving frequency (Y/N; PINT only supports N)",
+            ),
+            "",
+        )
+        self.add_param_from_top(
+            floatParameter(
+                name="DMDATA",
+                value=0.0,
+                units="",
+                description="Was the fit done using per-TOA DM information?",
+            ),
+            "",
+        )
+        self.add_param_from_top(
+            floatParameter(
+                name="NTOA",
+                value=0.0,
+                units="",
+                description="Number of TOAs used in the fitting",
+            ),
+            "",
+        )
+        self.add_param_from_top(
+            floatParameter(
+                name="CHI2",
+                value=0.0,
+                units="",
+                description="Chi-squared value obtained during fitting",
+            ),
+            "",
+        )
 
         for cp in components:
             self.add_component(cp, validate=False)
@@ -246,6 +296,14 @@ def validate(self):
 
         The checks include required parameters and parameter values.
         """
+        if self.DILATEFREQ.value:
+            warn("PINT does not support 'DILATEFREQ Y'")
+        if self.TIMEEPH.value not in [None, "FB90"]:
+            warn("PINT only supports 'TIMEEPH FB90'")
+        if self.T2CMETHOD.value not in [None, "IAU2000B"]:  # FIXME: really?
+            warn("PINT only supports 'T2CMETHOD IAU2000B'")
+        if self.UNITS.value not in [None, "TDB"]:
+            raise ValueError("PINT only supports 'UNITS TDB'")
         for cp in self.components.values():
             cp.validate()
 
@@ -1875,6 +1933,7 @@ def as_parfile(
         last_order=["jump_delay"],
     ):
         """Represent the entire model as a parfile string."""
+        self.validate()
         result_begin = ""
         result_end = ""
         result_middle = ""

tests/test_fitter_error_checking.py

@@ -4,6 +4,7 @@
 import numpy as np
 import astropy.units as u
 import pytest
+import re
 
 import pint.fitter
 from pint.models import get_model
@@ -17,6 +18,7 @@
 ELONG 0 0
 PEPOCH 57000
 DM 10 0
+SOLARN0 0
 """
 
 
@@ -130,3 +132,29 @@ def test_jump_everything_wideband():
         fitter.fit_toas(threshold=1e-14)
     for p in fitter.model.free_params:
         assert not np.isnan(fitter.model[p].value)
+
+
+@pytest.mark.parametrize("Fitter", [pint.fitter.WLSFitter, pint.fitter.GLSFitter])
+def test_update_model(Fitter):
+    model = get_model(io.StringIO("\n".join([par_base, "JUMP TEL barycenter 0"])))
+    model.INFO.value = "-f"
+    model.ECL.value = "IERS2010"
+    model.TIMEEPH.value = "FB90"
+    model.T2CMETHOD.value = "IERS2000B"
+    toas = make_fake_toas(58000, 59000, 10, model, obs="barycenter", freq=np.inf)
+    fitter = Fitter(toas, model)
+    fitter.fit_toas()
+    par_out = fitter.model.as_parfile()
+    assert re.search(r"CLOCK *TT\(TAI\)", par_out)
+    assert re.search(r"TIMEEPH *FB90", par_out)
+    assert re.search(r"T2CMETHOD *IERS2000B", par_out)
+    assert re.search(r"NE_SW *0.0", par_out)
+    assert re.search(r"ECL *IERS2010", par_out)
+    assert re.search(r"DILATEFREQ *N", par_out)
+    assert re.search(r"INFO *-f", par_out)
+    assert re.search(r"NTOA *10.0", par_out)
+    assert re.search(r"CHI2 *\d+.\d+", par_out)
+    assert re.search(r"EPHEM *DE421", par_out)
+    assert re.search(r"DMDATA *0.0", par_out)
+    assert re.search(r"START *58000.0", par_out)
+    assert re.search(r"FINISH *59000.0", par_out)