From 5c880f2295646f8299a3f49c7757e4b7f10f62cf Mon Sep 17 00:00:00 2001
From: Scott Ransom <sransom@nrao.edu>
Date: Fri, 8 Jul 2022 23:43:24 -0400
Subject: [PATCH 1/8] Improvements to designmatrix and fast random models

---
 src/pint/models/timing_model.py | 118 ++++++++++----
 src/pint/pintk/plk.py           |   6 +-
 src/pint/pintk/pulsar.py        |   4 +-
 src/pint/residuals.py           |  65 ++------
 src/pint/simulation.py          | 268 ++++++++++++++++++++++----------
 5 files changed, 291 insertions(+), 170 deletions(-)

diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py
index 69640e69a..e178a5b4c 100644
--- a/src/pint/models/timing_model.py
+++ b/src/pint/models/timing_model.py
@@ -60,9 +60,12 @@
 from pint.toa import TOAs
 from pint.utils import (
     PrefixError,
+    interesting_lines,
+    lines_of,
     split_prefixed_name,
     open_or_use,
     colorize,
+    taylor_horner
 )
 
 
@@ -1531,6 +1534,57 @@ def delete_jump_and_flags(self, toa_table, jump_num):
             return
         self.components["PhaseJump"].setup()
 
+    def get_spin_freq(self, times, calctype="modelF0"):
+        """Return barycentric pulsar rotational frequency in Hz at specific times
+
+        Parameters
+        ----------
+        times : float or long double MJD (can be array), astropy.Time object, TOAs
+            The times (barycentric, inf freq) at which to compute the spin frequency
+        calctype : {'modelF0', 'numerical', 'taylor'}
+            Type of calculation:
+                "modelF0" means a constant ``F0`` from the timing model
+                "numerical" uses the numerical derivative ``d_phase_d_toa()``
+                    This option requires times to be a TOAs object
+                "taylor" uses a Taylor series to compute the frequencies
+
+        Returns
+        -------
+        freq : astropy.units.Quantity (can be array)
+            The spin frequency in Hz at each barycentric time
+
+        Note
+        ----
+        "numerical" is significantly slower, but is much more exact, and
+        """
+        calc = calctype.lower()
+        assert calc in ["modelf0", "taylor", "numerical"]
+
+        if calc == "modelf0":
+            return self.F0.quantity
+
+        # ToDo: How does this handle Glitch or Piecewise models?
+        if calc == "numerical":
+            assert isinstance(times, TOAs)
+            return self.d_phase_d_toa(times)
+
+        # Handle various types of input times for "taylor" calc
+        if isinstance(times, TOAs):
+            # If we pass the TOA table, then barycenter the TOAs
+            bts = self.get_barycentric_toas(times)  # have units of days
+        elif isinstance(times, time.Time):
+            if times.scale == "tdb":
+                bts = np.asarray(times.mjd_long) << u.d
+            else:
+                raise ValueError("times as Time instance in needs scale=='tdb'")
+        elif isinstance(times, MJDParameter):
+            bts = np.asarray(times.value) << u.d  # .value is always a MJD long double
+        else:
+            bts = np.asarray(times) << u.d
+        dts = bts - (self.PEPOCH.value << u.d)
+        fterms = self.get_spin_terms()
+        return taylor_horner(dts, fterms).to(u.Hz)
+
     def get_barycentric_toas(self, toas, cutoff_component=""):
         """Conveniently calculate the barycentric TOAs.
 
@@ -1749,14 +1803,16 @@ def d_dm_d_param(self, data, param):
             )
         return result
 
-    def designmatrix(self, toas, acc_delay=None, incfrozen=False, incoffset=True):
+    def designmatrix(self, toas, incfrozen=False, incoffset=True, calctype="taylor"):
         """Return the design matrix.
 
-        The design matrix is the matrix with columns of ``d_phase_d_param/F0``
-        or ``d_toa_d_param``; it is used in fitting and calculating parameter
-        covariances.
+        The design matrix is the matrix with columns of
+        ``d_phase_d_param/F0`` or ``d_toa_d_param`` (or possibly
+        ``d_phase_d_param`` if calctype==``phase``); it is used in fitting
+        and calculating parameter covariances.
 
-        The value of ``F0`` used here is the parameter value in the model.
+        If calctype is not ``phase``, then it computes the local spin freq
+        by default via the ``taylor`` method.
 
         The order of parameters that are included is that returned by
         ``self.params``.
@@ -1765,29 +1821,36 @@ def designmatrix(self, toas, acc_delay=None, incfrozen=False, incoffset=True):
         ----------
         toas : pint.toa.TOAs
             The TOAs at which to compute the design matrix.
-        acc_delay
-            ???
         incfrozen : bool
             Whether to include frozen parameters in the design matrix
         incoffset : bool
             Whether to include the constant offset in the design matrix
+        calctype : {'modelF0', 'numerical', 'taylor'}
+            Type of calculation (for ``.get_spin_freq()``):
+                "modelF0" means a constant ``F0`` from the timing model
+                "numerical" uses the numerical derivative ``d_phase_d_toa()``
+                    This option requires times to be a TOAs object
+                "taylor" uses a Taylor series to compute the frequencies
+                "phase" return the matrix in terms of phase, not time
 
         Returns
         -------
         M : array
-            The design matrix, with shape (len(toas), len(self.free_params)+1)
-        names : list of str
-            The names of parameters in the corresponding parts of the design matrix
+            The design matrix, with shape (len(toas),
+            len(self.free_params)+incoffset)
+        names : list of strings
+            The names of parameters in the corresponding parts of the
+            design matrix
         units : astropy.units.Unit
             The units of the corresponding parts of the design matrix
 
         Note
         ----
-        Here we have negative sign here. Since in pulsar timing
-        the residuals are calculated as (Phase - int(Phase)), which is different
-        from the conventional definition of least square definition (Data - model)
-        We decide to add minus sign here in the design matrix, so the fitter
-        keeps the conventional way.
+        We return the negative of the ``d_phase_d_param()`` since in
+        pulsar timing the residuals are calculated as (Phase -
+        int(Phase)), which is different from the conventional definition
+        of least square definition (Data - model).  Including the negative
+        here allows the fitter to keep the standard convention.
         """
 
         noise_params = self.get_params_of_component_type("NoiseComponent")
@@ -1805,26 +1868,27 @@ def designmatrix(self, toas, acc_delay=None, incfrozen=False, incoffset=True):
             if (incfrozen or not getattr(self, par).frozen) and par not in noise_params
         ]
 
-        F0 = self.F0.quantity  # 1/sec
-        ntoas = len(toas)
-        nparams = len(params)
+        # This is either ones if we want the derivs in phase rather than
+        # time, or it is the spin freq calculated by the chosen method
+        norm = (
+            self.get_spin_freq(toas, calctype)
+            if calctype != "phase"
+            else (np.ones(toas.ntoas) << u.dimensionless_unscaled)
+        )
         delay = self.delay(toas)
         units = []
-        # Apply all delays ?
-        # tt = toas['tdbld']
-        # for df in self.delay_funcs:
-        #    tt -= df(toas)
 
-        M = np.zeros((ntoas, nparams))
+        # SMR: Wouldn't it be better to have this in thetransposed shape?
+        M = np.zeros((toas.ntoas, len(params)), dtype=float)
         for ii, param in enumerate(params):
             if param == "Offset":
-                M[:, ii] = 1.0 / F0.value
-                units.append(u.s / u.s)
+                M[:, ii] = 1.0 / norm.value
+                units.append(u.Unit("") / norm.unit)
             else:
                 q = -self.d_phase_d_param(toas, delay, param)
                 the_unit = u.Unit("") / getattr(self, param).units
-                M[:, ii] = q.to_value(the_unit) / F0.value
-                units.append(the_unit / F0.unit)
+                M[:, ii] = q.to_value(the_unit) / norm.value
+                units.append(the_unit / norm.unit)
         return M, params, units
 
     def compare(
diff --git a/src/pint/pintk/plk.py b/src/pint/pintk/plk.py
index 06aee245e..43e0bd24f 100644
--- a/src/pint/pintk/plk.py
+++ b/src/pint/pintk/plk.py
@@ -1070,11 +1070,11 @@ def plotResiduals(self, keepAxes=False):
                     else self.psr.postfit_resids
                 )
                 if self.y_unit == u.us:
-                    f0 = r.get_PSR_freq().to(u.MHz).value
+                    f0 = r.model.F0.quantity.to(u.MHz).value
                 elif self.y_unit == u.ms:
-                    f0 = r.get_PSR_freq().to(u.kHz).value
+                    f0 = r.model.F0.quantity.to(u.kHz).value
                 else:
-                    f0 = r.get_PSR_freq().to(u.Hz).value
+                    f0 = r.model.F0.quantity.to(u.Hz).value
                 self.plkAx2x.set_visible(True)
                 self.plkAx2x.set_ylabel(plotlabels[self.yid][1])
                 self.plkAx2x.set_ylim(ymin * f0, ymax * f0)
diff --git a/src/pint/pintk/pulsar.py b/src/pint/pintk/pulsar.py
index 58e33e323..8584d0b87 100644
--- a/src/pint/pintk/pulsar.py
+++ b/src/pint/pintk/pulsar.py
@@ -632,7 +632,7 @@ def random_models(self, selected):
                 Ntoas,
                 self.postfit_model,
                 obs="coe",
-                freq=1 * u.THz,  # effectively infinite frequency
+                freq=np.median(self.all_toas.get_freqs()),
                 include_bipm=sim_sel.clock_corr_info["include_bipm"],
                 include_gps=sim_sel.clock_corr_info["include_gps"],
             )
@@ -653,6 +653,8 @@ def random_models(self, selected):
             Nmodels=15,
             keep_models=False,
             return_time=True,
+            params="all",
+            fast_method=True,
         )
 
         # Get a selection array for the fake TOAs that covers the fit TOAs (plus extra)
diff --git a/src/pint/residuals.py b/src/pint/residuals.py
index a80b323ab..41d8a0759 100644
--- a/src/pint/residuals.py
+++ b/src/pint/residuals.py
@@ -243,54 +243,6 @@ def rms_weighted(self):
         wmean, werr, wsdev = weighted_mean(self.time_resids, w, sdev=True)
         return wsdev.to(u.us)
 
-    def get_PSR_freq(self, calctype="modelF0"):
-        """Return pulsar rotational frequency in Hz.
-
-        Parameters
-        ----------
-        calctype : {'modelF0', 'numerical', 'taylor'}
-            Type of calculation.  If `calctype` == "modelF0", then simply the ``F0``
-            parameter from the model.
-            If `calctype` == "numerical", then try a numerical derivative
-            If `calctype` == "taylor", evaluate the frequency with a Taylor series
-
-        Returns
-        -------
-        freq : astropy.units.Quantity
-            Either the single ``F0`` in the model or the spin frequency at the moment of each TOA.
-        """
-        assert calctype.lower() in ["modelf0", "taylor", "numerical"]
-        if calctype.lower() == "modelf0":
-            # TODO this function will be re-write and move to timing model soon.
-            # The following is a temproary patch.
-            if "Spindown" in self.model.components:
-                F0 = self.model.F0.quantity
-            elif "P0" in self.model.params:
-                F0 = 1.0 / self.model.P0.quantity
-            else:
-                raise AttributeError(
-                    "No pulsar spin parameter(e.g., 'F0'," " 'P0') found."
-                )
-            return F0.to(u.Hz)
-        elif calctype.lower() == "taylor":
-            # see Spindown.spindown_phase
-            dt = self.model.get_dt(self.toas, 0)
-            # if the model is defined through F0, F1, ...
-            if "F0" in self.model.params:
-                fterms = [0.0 * u.dimensionless_unscaled] + self.model.get_spin_terms()
-
-            # otherwise assume P0, PDOT
-            else:
-                F0 = 1.0 / self.model.P0.quantity
-                if "PDOT" in self.model.params:
-                    F1 = -self.model.PDOT.quantity / self.model.P0.quantity**2
-                else:
-                    F1 = 0 * u.Hz / u.s
-                fterms = [0.0 * u.dimensionless_unscaled, F0, F1]
-            return taylor_horner_deriv(dt, fterms, deriv_order=1).to(u.Hz)
-        elif calctype.lower() == "numerical":
-            return self.model.d_phase_d_toa(self.toas)
-
     def calc_phase_resids(self):
         """Compute timing model residuals in pulse phase."""
 
@@ -369,11 +321,12 @@ def calc_time_resids(self, calctype="taylor"):
 
         Parameters
         ----------
-        calctype : {'taylor', 'modelF0', 'numerical'}
-            Type of calculation.  If `calctype` == "modelF0", then simply the ``F0``
-            parameter from the model.
-            If `calctype` == "numerical", then try a numerical derivative
-            If `calctype` == "taylor", evaluate the frequency with a Taylor series
+        calctype : {'modelF0', 'numerical', 'taylor'}
+            Type of calculation:
+                "modelF0" means a constant ``F0`` from the timing model
+                "numerical" uses the numerical derivative ``d_phase_d_toa()``
+                    This option requires times to be a TOAs object
+                "taylor" uses a Taylor series to compute the frequencies
 
         Returns
         -------
@@ -381,12 +334,14 @@ def calc_time_resids(self, calctype="taylor"):
 
         See Also
         --------
-        :meth:`pint.residuals.Residuals.get_PSR_freq`
+        :meth:`pint.timing_model.TimingModel.get_spin_freq`
         """
         assert calctype.lower() in ["modelf0", "taylor", "numerical"]
         if self.phase_resids is None:
             self.phase_resids = self.calc_phase_resids()
-        return (self.phase_resids / self.get_PSR_freq(calctype=calctype)).to(u.s)
+        return (
+            self.phase_resids / self.model.get_spin_freq(self.toas, calctype=calctype)
+        ).to(u.s)
 
     def calc_chi2(self, full_cov=False):
         """Return the weighted chi-squared for the model and toas.
diff --git a/src/pint/simulation.py b/src/pint/simulation.py
index 65bba9be4..45c0e8ef6 100644
--- a/src/pint/simulation.py
+++ b/src/pint/simulation.py
@@ -8,9 +8,10 @@
 from loguru import logger as log
 from astropy import time
 
-import pint.residuals
+import pint.residuals as pr
 import pint.toa
-from pint.observatory import bipm_default, get_observatory
+import pint.models
+from pint.observatory import Observatory, bipm_default, get_observatory
 
 __all__ = [
     "zero_residuals",
@@ -69,7 +70,7 @@ def zero_residuals(ts, model, maxiter=10, tolerance=None):
         else:
             tolerance = 5 * u.us
     for i in range(maxiter):
-        r = pint.residuals.Residuals(ts, model, track_mode="use_pulse_numbers")
+        r = pr.Residuals(ts, model, track_mode="use_pulse_numbers")
         resids = r.calc_time_resids(calctype="taylor")
         if maxresid is not None and (np.abs(resids).max() > maxresid):
             log.warning(
@@ -385,76 +386,36 @@ def make_fake_toas_fromtim(timfile, model, add_noise=False, name="fake"):
     return make_fake_toas(input_ts, model=model, add_noise=add_noise, name=name)
 
 
-def calculate_random_models(
-    fitter, toas, Nmodels=100, keep_models=True, return_time=False, params="all"
-):
+def make_random_models(fitter, Nmodels=100, params="all"):
     """
-    Calculates random models based on the covariance matrix of the `fitter` object.
-
-    returns the new phase differences compared to the original model
-    optionally returns all of the random models
+    Create random models based on the covariance matrix of fitter
 
     Parameters
     ----------
-    fitter: pint.fitter.Fitter
-        current fitter object containing a model and parameter covariance matrix
-    toas: pint.toa.TOAs
-        TOAs to calculate models
-    Nmodels: int, optional
-        number of random models to calculate
-    keep_models: bool, optional
-        whether to keep and return the individual random models (slower)
-    params: list, optional
-        if specified, selects only those parameters to vary.  Default ('all') is to use all parameters other than Offset
+    fitter : pint.fitter.Fitter
+        current post-fit fitter object instance containing a model
+    Nmodels : int, default=100
+        number of random models to create
+    params : list, default="all"
+        if specified, selects only those parameters to vary.  Default
+        ('all') is to use all parameters other than Offset
 
     Returns
     -------
-    dphase : np.ndarray
-        phase difference with respect to input model, size is [Nmodels, len(toas)]
-    random_models : list, optional
+    random_models : list
         list of random models (each is a :class:`pint.models.timing_model.TimingModel`)
 
-    Example
-    -------
-    >>> from pint.models import get_model_and_toas
-    >>> from pint import fitter, toa
-    >>> import pint.simulation
-    >>> import io
-    >>>
-    >>> # the locations of these may vary
-    >>> timfile = "tests/datafile/NGC6440E.tim"
-    >>> parfile = "tests/datafile/NGC6440E.par"
-    >>> m, t = get_model_and_toas(parfile, timfile)
-    >>> # fit the model to the data
-    >>> f = fitter.WLSFitter(toas=t, model=m)
-    >>> f.fit_toas()
-    >>>
-    >>> # make fake TOAs starting at the end of the
-    >>> # current data and going out 100 days
-    >>> tnew = simulation.make_fake_toas_uniform(t.get_mjds().max().value,
-    >>>                           t.get_mjds().max().value+100, 50, model=f.model)
-    >>> # now make random models
-    >>> dphase, mrand = pint.simulation.calculate_random_models(f, tnew, Nmodels=100)
-
-
     Note
     ----
-    To calculate new TOAs, you can use :func:`~pint.simulation.make_fake_toas`
-
-    or similar
+    This function is used by :func:`~pint.simulation.calculate_random_models`
     """
-    Nmjd = len(toas)
-    phases_i = np.zeros((Nmodels, Nmjd))
-    phases_f = np.zeros((Nmodels, Nmjd))
-    freqs = np.zeros((Nmodels, Nmjd), dtype=np.float128) * u.Hz
-
     cov_matrix = fitter.parameter_covariance_matrix
-    # this is a list of the parameter names in the order they appear in the coviarance matrix
+    # this is a list of the parameter names in the order they appear in the covariance matrix
     param_names = cov_matrix.get_label_names(axis=0)
     # this is a dictionary with the parameter values, but it might not be in the same order
     # and it leaves out the Offset parameter
     param_values = fitter.model.get_params_dict("free", "value")
-    mean_vector = np.array([param_values[x] for x in param_names if not x == "Offset"])
+    mean_vector = np.array([param_values[x] for x in param_names if x != "Offset"])
     if params == "all":
         # remove the first column and row (absolute phase)
         if param_names[0] == "Offset":
@@ -478,36 +439,175 @@ def calculate_random_models(
             mean_vector = mean_vector[index]
         param_names = cov_matrix.get_label_names(axis=0)
 
-    f_rand = deepcopy(fitter)
-
     # scale by fac
     mean_vector = mean_vector * fac
     scaled_cov_matrix = ((cov_matrix.matrix * fac).T * fac).T
     random_models = []
-    for imodel in range(Nmodels):
+    for _ in range(Nmodels):
+        f_rand = deepcopy(fitter)
         # create a set of randomized parameters based on mean vector and covariance matrix
-        rparams_num = np.random.multivariate_normal(mean_vector, scaled_cov_matrix)
-        # scale params back to real units
-        for j in range(len(mean_vector)):
-            rparams_num[j] /= fac[j]
-        rparams = OrderedDict(zip(param_names, rparams_num))
-        f_rand.set_params(rparams)
-        phase = f_rand.model.phase(toas, abs_phase=True)
-        phases_i[imodel] = phase.int
-        phases_f[imodel] = phase.frac
-        r = pint.residuals.Residuals(toas, f_rand.model)
-        freqs[imodel] = r.get_PSR_freq(calctype="taylor")
-        if keep_models:
-            random_models.append(f_rand.model)
-            f_rand = deepcopy(fitter)
-    phases = phases_i + phases_f
-    phases0 = fitter.model.phase(toas, abs_phase=True)
-    dphase = phases - (phases0.int + phases0.frac)
-
-    if return_time:
-        dphase /= freqs
-
-    if keep_models:
-        return dphase, random_models
-    else:
-        return dphase
+        # dividing by fac brings us back to real units
+        rparams = np.random.multivariate_normal(mean_vector, scaled_cov_matrix) / fac
+        f_rand.set_params(dict(zip(param_names, rparams)))
+        random_models.append(f_rand.model)
+    return random_models
+
+
+def compute_random_model_resids_exact(fitter, models, toas, return_time=False):
+    """
+    Calculate residuals from random models exactly
+
+    Parameters
+    ----------
+    fitter : pint.fitter.Fitter
+        current post-fit fitter object instance containing the reference model
+    models : list, pint.models.timing_model.TimingModel
+        list of the random models to compute residuals
+    toas : pint.toa.TOAs
+        TOAs to calculate residuals from
+    return_time : bool, default=False
+        return the residuals in time rather than in phase
+
+    Returns
+    -------
+    dphase : np.ndarray
+        phase difference with respect to reference model, size is [len(models), len(toas)]
+
+    Note
+    ----
+    This function is used by :func:`~pint.simulation.calculate_random_models`
+    """
+    resids = np.zeros((len(models), toas.ntoas), dtype=np.float)
+    # These are the reference residuals from the reference model
+    r0 = pr.Residuals(toas, fitter.model, subtract_mean=False)
+    for ii, model in enumerate(models):
+        rn = pr.Residuals(toas, model, subtract_mean=False)
+        resids[ii] = (
+            rn.time_resids - r0.time_resids
+            if return_time
+            else rn.phase_resids - r0.phase_resids
+        )
+    return resids << u.s if return_time else resids
+
+
+def compute_random_model_resids_fast(fitter, models, toas, return_time=False):
+    """
+    Calculate residuals from random models using the design matrix
+
+    Parameters
+    ----------
+    fitter : pint.fitter.Fitter
+        current post-fit fitter object instance containing the reference model
+    models : list, pint.models.timing_model.TimingModel
+        list of the random models to compute residuals
+    toas : pint.toa.TOAs
+        TOAs to calculate residuals from
+    return_time : bool, default=False
+        return the residuals in time rather than in phase
+
+    Returns
+    -------
+    dphase : np.ndarray
+        phase difference with respect to reference model, size is [len(models), len(toas)]
+
+    Note
+    ----
+    This function is used by :func:`~pint.simulation.calculate_random_models`.  It should be
+    a factor of about three faster than the exact version, and the majority of the time it
+    takes is calculating the design matrix.  The residual computation via np.dot() is
+    incredibly fast.
+    """
+    # Compute the design matrix for the toas using the reference model
+    m = fitter.model
+    calctype = "taylor" if return_time else "phase"
+    dmat, labels, _ = m.designmatrix(
+        toas, incfrozen=False, incoffset=False, calctype=calctype
+    )
+    rs = np.zeros((len(models), toas.ntoas), dtype=np.float)
+    dparam = np.zeros(len(labels), dtype=np.float128)
+    for ii, rm in enumerate(models):
+        # Do the differencing between the random models and the reference model
+        for jj, pp in enumerate(labels):
+            if type(getattr(m, pp)) == pint.models.parameter.MJDParameter:
+                dparam[jj] = (getattr(m, pp).value - getattr(rm, pp).value) * u.d
+            else:
+                dparam[jj] = getattr(m, pp).quantity - getattr(rm, pp).quantity
+        # Compute residuals via a dot product of the param diffs with the design matrix
+        rs[ii] = np.dot(dparam, dmat.T)
+    return rs << u.s if return_time else rs
+
+
+def calculate_random_models(
+    fitter,
+    toas,
+    Nmodels=100,
+    keep_models=True,
+    return_time=False,
+    params="all",
+    fast_method=True,
+):
+    """
+    Calculates random model residuals based on the parameter covariance matrix
+
+    Parameters
+    ----------
+    fitter : pint.fitter.Fitter
+        current post-fit fitter object instance containing a model
+    toas : pint.toa.TOAs
+        TOAs to calculate models
+    Nmodels : int, default=100
+        number of random models to calculate
+    keep_models : bool, default=True
+        whether to keep and return the individual random models (slower)
+    return_time : bool, default=False
+        return the residuals in time rather than in phase
+    params : list, default="all"
+        if specified, selects only those parameters to vary.  Default
+        ('all') is to use all parameters other than Offset
+    fast_method : bool, default=True
+        Compute the residuals using the design matrix and perturbations
+        rather than full residual calculation and subtraction from the
+        reference model (typically 3-5x faster)
+
+    Returns
+    -------
+    dphase : np.ndarray
+        phase difference with respect to input model, size is [Nmodels, len(toas)]
+    random_models : list, optional
+        list of random models (each is a :class:`pint.models.timing_model.TimingModel`)
+
+    Example
+    -------
+    >>> from pint.models import get_model_and_toas
+    >>> from pint import fitter, toa
+    >>> import pint.simulation
+    >>> import io
+    >>>
+    >>> # the locations of these may vary
+    >>> timfile = "tests/datafile/NGC6440E.tim"
+    >>> parfile = "tests/datafile/NGC6440E.par"
+    >>> m, t = get_model_and_toas(parfile, timfile)
+    >>> # fit the model to the data
+    >>> f = fitter.WLSFitter(toas=t, model=m)
+    >>> f.fit_toas()
+    >>>
+    >>> # make fake TOAs starting at the end of the
+    >>> # current data and going out 100 days
+    >>> tnew = simulation.make_fake_toas_uniform(t.get_mjds().max().value,
+    >>>                           t.get_mjds().max().value+100, 50, model=f.model)
+    >>> # now make random models
+    >>> dphase, mrand = pint.simulation.calculate_random_models(f, tnew, Nmodels=100)
+
+    Note
+    ----
+    To calculate new TOAs, you can use
+    :func:`~pint.simulation.make_fake_toas` or similar
+    """
+    models = make_random_models(fitter, Nmodels, params)
+    residuals_fn = (
+        compute_random_model_resids_fast
+        if fast_method
+        else compute_random_model_resids_exact
+    )
+    resids = residuals_fn(fitter, models, toas, return_time)
+    return (resids, models) if keep_models else resids

From 9573efcf5ccefbba017559890032a4cd678bb34c Mon Sep 17 00:00:00 2001
From: Scott Ransom <sransom@nrao.edu>
Date: Sat, 9 Jul 2022 10:27:18 -0400
Subject: [PATCH 2/8] Changed PEPOCH_P0 to PEPOCH and added get_spin_terms()
 method

---
 docs/examples/How_to_build_a_timing_model_component.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/docs/examples/How_to_build_a_timing_model_component.py b/docs/examples/How_to_build_a_timing_model_component.py
index 0dfe35264..2dbb5657d 100644
--- a/docs/examples/How_to_build_a_timing_model_component.py
+++ b/docs/examples/How_to_build_a_timing_model_component.py
@@ -123,7 +123,7 @@ def __init__(self):
         # Add reference epoch time.
         self.add_param(
             p.MJDParameter(
-                name="PEPOCH_P0",
+                name="PEPOCH",
                 description="Reference epoch for spin-down",
                 time_scale="tdb",
             )
@@ -188,7 +188,11 @@ def d_F1_d_P1(self):
     def get_dt(self, toas, delay):
         """dt from the toas to the reference time."""
         # toas.table['tdbld'] stores the tdb time in longdouble.
-        return (toas.table["tdbld"] - self.PEPOCH_P0.value) * u.day - delay
+        return (toas.table["tdbld"] - self.PEPOCH.value) * u.day - delay
+
+    def get_spin_terms(self):
+        """Return a list of the spin term values in the model: [F0, F1, ..., FN]."""
+        return [self.F0.quantity, self.F1.quantity]
 
     # Defining the phase function, which is added to the self.phase_funcs_component
     def spindown_phase_period(self, toas, delay):
@@ -222,7 +226,7 @@ def d_phase_d_P1(self, toas, param, delay):
              DECJ      -20:21:29.0  1 0.4
              P0        0.016264003404474613 1 0
              P1        3.123955D-19 1 0
-             PEPOCH_P0     53750.000000
+             PEPOCH        53750.000000
              POSEPOCH      53750.000000
              DM              223.9  1 0.3
              SOLARN0               0.00

From 29e7f6809f688a21696e026726258467ed8c3940 Mon Sep 17 00:00:00 2001
From: Scott Ransom <sransom@nrao.edu>
Date: Sun, 10 Jul 2022 13:58:44 -0400
Subject: [PATCH 3/8] only copy ref model rather than full fitter when making
 random models

---
 src/pint/simulation.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/pint/simulation.py b/src/pint/simulation.py
index 45c0e8ef6..34fec0150 100644
--- a/src/pint/simulation.py
+++ b/src/pint/simulation.py
@@ -444,12 +444,12 @@ def make_random_models(fitter, Nmodels=100, params="all"):
     scaled_cov_matrix = ((cov_matrix.matrix * fac).T * fac).T
     random_models = []
     for _ in range(Nmodels):
-        f_rand = deepcopy(fitter)
+        m_rand = deepcopy(fitter.model)
         # create a set of randomized parameters based on mean vector and covariance matrix
         # dividing by fac brings us back to real units
         rparams = np.random.multivariate_normal(mean_vector, scaled_cov_matrix) / fac
-        f_rand.set_params(dict(zip(param_names, rparams)))
-        random_models.append(f_rand.model)
+        m_rand.set_param_values(dict(zip(param_names, rparams)))
+        random_models.append(m_rand)
     return random_models
 
 

From 2db9de110abb58edc5e1a1224c729903b46f0d14 Mon Sep 17 00:00:00 2001
From: Scott Ransom <sransom@nrao.edu>
Date: Sun, 10 Jul 2022 14:12:43 -0400
Subject: [PATCH 4/8] Fixed docstring first lines

---
 src/pint/simulation.py | 12 ++++--------
 1 file changed, 4 insertions(+), 8 deletions(-)

diff --git a/src/pint/simulation.py b/src/pint/simulation.py
index 34fec0150..6a1d6255b 100644
--- a/src/pint/simulation.py
+++ b/src/pint/simulation.py
@@ -387,8 +387,7 @@ def make_fake_toas_fromtim(timfile, model, add_noise=False, name="fake"):
 
 
 def make_random_models(fitter, Nmodels=100, params="all"):
-    """
-    Create random models based on the covariance matrix of fitter
+    """Create random models based on the covariance matrix of fitter.
 
     Parameters
     ----------
@@ -454,8 +453,7 @@ def make_random_models(fitter, Nmodels=100, params="all"):
 
 
 def compute_random_model_resids_exact(fitter, models, toas, return_time=False):
-    """
-    Calculate residuals from random models exactly
+    """Calculate residuals from random models exactly.
 
     Parameters
     ----------
@@ -491,8 +489,7 @@ def compute_random_model_resids_exact(fitter, models, toas, return_time=False):
 
 
 def compute_random_model_resids_fast(fitter, models, toas, return_time=False):
-    """
-    Calculate residuals from random models using the design matrix
+    """Calculate residuals from random models using the design matrix.
 
     Parameters
     ----------
@@ -546,8 +543,7 @@ def calculate_random_models(
     params="all",
     fast_method=True,
 ):
-    """
-    Calculates random model residuals based on the parameter covariance matrix
+    """Calculates random model residuals based on the parameter covariance matrix.
 
     Parameters
     ----------

From 110425fd7ccf5a0c68ee4ce92fddf23cc5caafc6 Mon Sep 17 00:00:00 2001
From: Scott Ransom <sransom@nrao.edu>
Date: Sun, 10 Jul 2022 14:28:54 -0400
Subject: [PATCH 5/8] Specify Residuals load directly

---
 src/pint/simulation.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/pint/simulation.py b/src/pint/simulation.py
index 6a1d6255b..528b77865 100644
--- a/src/pint/simulation.py
+++ b/src/pint/simulation.py
@@ -8,9 +8,9 @@
 from loguru import logger as log
 from astropy import time
 
-import pint.residuals as pr
 import pint.toa
 import pint.models
+from pint.residuals import Residuals
 from pint.observatory import Observatory, bipm_default, get_observatory
 
 __all__ = [
@@ -70,7 +70,7 @@ def zero_residuals(ts, model, maxiter=10, tolerance=None):
         else:
             tolerance = 5 * u.us
     for i in range(maxiter):
-        r = pr.Residuals(ts, model, track_mode="use_pulse_numbers")
+        r = Residuals(ts, model, track_mode="use_pulse_numbers")
         resids = r.calc_time_resids(calctype="taylor")
         if maxresid is not None and (np.abs(resids).max() > maxresid):
             log.warning(
@@ -477,9 +477,9 @@ def compute_random_model_resids_exact(fitter, models, toas, return_time=False):
     """
     resids = np.zeros((len(models), toas.ntoas), dtype=np.float)
     # These are the reference residuals from the reference model
-    r0 = pr.Residuals(toas, fitter.model, subtract_mean=False)
+    r0 = Residuals(toas, fitter.model, subtract_mean=False)
     for ii, model in enumerate(models):
-        rn = pr.Residuals(toas, model, subtract_mean=False)
+        rn = Residuals(toas, model, subtract_mean=False)
         resids[ii] = (
             rn.time_resids - r0.time_resids
             if return_time

From 145abb3ae7113c044c6ed59f5ec5da5fa3b597d6 Mon Sep 17 00:00:00 2001
From: Scott Ransom <sransom@nrao.edu>
Date: Sun, 10 Jul 2022 14:41:26 -0400
Subject: [PATCH 6/8] Replace two asserts with exception raises, if needed

---
 src/pint/models/timing_model.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py
index e178a5b4c..1652a4d9d 100644
--- a/src/pint/models/timing_model.py
+++ b/src/pint/models/timing_model.py
@@ -1558,14 +1558,18 @@ def get_spin_freq(self, times, calctype="modelF0"):
         "numerical" is significantly slower, but is much more exact, and
         """
         calc = calctype.lower()
-        assert calc in ["modelf0", "taylor", "numerical"]
+        if calc not in ["modelf0", "taylor", "numerical"]:
+            raise ValueError(
+                "calctype must be one of ['modelf0', 'taylor'', 'numerical']"
+            )
 
         if calc == "modelf0":
             return self.F0.quantity
 
         # ToDo: How does this handle Glitch or Piecewise models?
         if calc == "numerical":
-            assert isinstance(times, TOAs)
+            if not isinstance(times, TOAs):
+                raise TypeError("times must be TOAs when using calctype='numerical'")
             return self.d_phase_d_toa(times)
 
         # Handle various types of input times for "taylor" calc

From 58d77f8db756e3894007a4eee867363bfd8afffc Mon Sep 17 00:00:00 2001
From: Scott Ransom <sransom@nrao.edu>
Date: Tue, 15 Nov 2022 21:47:52 -0500
Subject: [PATCH 7/8] Blacked

---
 src/pint/models/timing_model.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py
index 1652a4d9d..5bd26e4c6 100644
--- a/src/pint/models/timing_model.py
+++ b/src/pint/models/timing_model.py
@@ -65,7 +65,7 @@
     split_prefixed_name,
     open_or_use,
     colorize,
-    taylor_horner
+    taylor_horner,
 )
 
 

From 00cf6f20c1ce5c9e50114192b324f2cef670fb8b Mon Sep 17 00:00:00 2001
From: Scott Ransom <sransom@nrao.edu>
Date: Tue, 15 Nov 2022 22:25:53 -0500
Subject: [PATCH 8/8] stylistic code cleanups recommended by Sorcery

---
 src/pint/models/timing_model.py | 210 +++++++++++++-------------------
 src/pint/simulation.py          |   7 +-
 2 files changed, 87 insertions(+), 130 deletions(-)

diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py
index 5bd26e4c6..025e1ffe6 100644
--- a/src/pint/models/timing_model.py
+++ b/src/pint/models/timing_model.py
@@ -60,8 +60,6 @@
 from pint.toa import TOAs
 from pint.utils import (
     PrefixError,
-    interesting_lines,
-    lines_of,
     split_prefixed_name,
     open_or_use,
     colorize,
@@ -88,22 +86,20 @@
 # errors in the par file.
 #
 # Comparisons with keywords in par file lines is done in a case insensitive way.
-ignore_params = set(
-    [
-        "TRES",
-        "TZRMJD",
-        "TZRFRQ",
-        "TZRSITE",
-        "NITS",
-        "IBOOT",
-        "CHI2R",
-        "MODE",
-        "PLANET_SHAPIRO2",
-        #    'NE_SW', 'NE_SW2',
-    ]
-)
-
-ignore_prefix = set(["DMXF1_", "DMXF2_", "DMXEP_"])  # DMXEP_ for now.
+ignore_params = {
+    "TRES",
+    "TZRMJD",
+    "TZRFRQ",
+    "TZRSITE",
+    "NITS",
+    "IBOOT",
+    "CHI2R",
+    "MODE",
+    "PLANET_SHAPIRO2",
+    # 'NE_SW', 'NE_SW2',
+}
+
+ignore_prefix = {"DMXF1_", "DMXF2_", "DMXEP_"}  # DMXEP_ for now.
 
 DEFAULT_ORDER = [
     "astrometry",
@@ -442,34 +438,29 @@ def params_ordered(self):
         used_cats = set()
         pstart = copy.copy(self.top_level_params)
         for cat in start_order:
-            if cat in compdict:
-                cp = compdict[cat]
-                for cpp in cp:
-                    pstart += cpp.params
-                used_cats.add(cat)
-            else:
+            if cat not in compdict:
                 continue
-
+            cp = compdict[cat]
+            for cpp in cp:
+                pstart += cpp.params
+            used_cats.add(cat)
         pend = []
         for cat in last_order:
-            if cat in compdict:
-                cp = compdict[cat]
-                for cpp in cp:
-                    pend += cpp.parms
-                used_cats.add(cat)
-            else:
+            if cat not in compdict:
                 continue
-
+            cp = compdict[cat]
+            for cpp in cp:
+                pend += cpp.parms
+            used_cats.add(cat)
         # Now collect any components that haven't already been included in the list
         pmid = []
         for cat in compdict:
             if cat in used_cats:
                 continue
-            else:
-                cp = compdict[cat]
-                for cpp in cp:
-                    pmid += cpp.params
-                used_cats.add(cat)
+            cp = compdict[cat]
+            for cpp in cp:
+                pmid += cpp.params
+            used_cats.add(cat)
 
         return pstart + pmid + pend
 
@@ -573,7 +564,7 @@ def get_params_of_component_type(self, component_type):
         -------
         list
         """
-        component_type_list_str = "{}_list".format(component_type)
+        component_type_list_str = f"{component_type}_list"
         if hasattr(self, component_type_list_str):
             component_type_list = getattr(self, component_type_list_str)
             return [
@@ -604,17 +595,14 @@ def set_param_uncertainties(self, fitp):
         """Set the model parameters to the value contained in the input dict."""
         for k, v in fitp.items():
             p = getattr(self, k)
-            if isinstance(v, u.Quantity):
-                p.uncertainty = v
-            else:
-                p.uncertainty = v * p.units
+            p.uncertainty = v if isinstance(v, u.Quantity) else v * p.units
 
     @property_exists
     def components(self):
         """All the components in a dictionary indexed by name."""
         comps = {}
         for ct in self.component_types:
-            for cp in getattr(self, ct + "_list"):
+            for cp in getattr(self, f"{ct}_list"):
                 comps[cp.__class__.__name__] = cp
         return comps
 
@@ -705,10 +693,7 @@ def orbital_phase(self, barytimes, anom="mean", radians=True):
             raise ValueError("anom='%s' is not a recognized type of anomaly" % anom)
         # Make sure all angles are between 0-2*pi
         anoms = np.remainder(anoms.value, 2 * np.pi)
-        if radians:  # return with radian units
-            return anoms * u.rad
-        else:  # return as unitless cycles from 0-1
-            return anoms / (2 * np.pi)
+        return anoms * u.rad if radians else anoms / (2 * np.pi)
 
     def conjunction(self, baryMJD):
         """Return the time(s) of the first superior conjunction(s) after baryMJD.
@@ -767,10 +752,7 @@ def funct(t):
                     break
             # Now use scipy to find the root
             scs.append(brentq(funct, ts[lb], ts[lb + 1]))
-        if len(scs) == 1:
-            return scs[0]  # Return a float
-        else:
-            return np.asarray(scs)  # otherwise return an array
+        return scs[0] if len(scs) == 1 else np.asarray(scs)
 
     @property_exists
     def dm_funcs(self):
@@ -868,11 +850,11 @@ def get_deriv_funcs(self, component_type, derivative_type=""):
         """Return dictionary of derivative functions."""
         # TODO, this function can be a more generical function collector.
         deriv_funcs = defaultdict(list)
-        if not derivative_type == "":
+        if derivative_type != "":
             derivative_type += "_"
-        for cp in getattr(self, component_type + "_list"):
+        for cp in getattr(self, f"{component_type}_list"):
             try:
-                df = getattr(cp, derivative_type + "deriv_funcs")
+                df = getattr(cp, f"{derivative_type}deriv_funcs")
             except AttributeError:
                 continue
             for k, v in df.items():
@@ -951,15 +933,14 @@ def map_component(self, component):
             if component not in list(comps.keys()):
                 raise AttributeError("No '%s' in the timing model." % component)
             comp = comps[component]
-        else:  # When component is an component instance.
-            if component not in list(comps.values()):
-                raise AttributeError(
-                    "No '%s' in the timing model." % component.__class__.__name__
-                )
-            else:
-                comp = component
+        elif component in list(comps.values()):
+            comp = component
+        else:
+            raise AttributeError(
+                "No '%s' in the timing model." % component.__class__.__name__
+            )
         comp_type = self.get_component_type(comp)
-        host_list = getattr(self, comp_type + "_list")
+        host_list = getattr(self, f"{comp_type}_list")
         order = host_list.index(comp)
         return comp, order, host_list, comp_type
 
@@ -978,9 +959,9 @@ def add_component(
             If true, add a duplicate component. Default is False.
         """
         comp_type = self.get_component_type(component)
+        cur_cps = []
         if comp_type in self.component_types:
-            comp_list = getattr(self, comp_type + "_list")
-            cur_cps = []
+            comp_list = getattr(self, f"{comp_type}_list")
             for cp in comp_list:
                 # If component order is not defined.
                 cp_order = (
@@ -1001,21 +982,19 @@ def add_component(
                     )
         else:
             self.component_types.append(comp_type)
-            cur_cps = []
-
         # link new component to TimingModel
         component._parent = self
 
-        # If the categore is not in the order list, it will be added to the end.
+        # If the category is not in the order list, it will be added to the end.
         if component.category not in order:
-            new_cp = tuple((len(order) + 1, component))
+            new_cp = len(order) + 1, component
         else:
-            new_cp = tuple((order.index(component.category), component))
+            new_cp = order.index(component.category), component
         # add new component
         cur_cps.append(new_cp)
         cur_cps.sort(key=lambda x: x[0])
         new_comp_list = [c[1] for c in cur_cps]
-        setattr(self, comp_type + "_list", new_comp_list)
+        setattr(self, f"{comp_type}_list", new_comp_list)
         # Set up components
         if setup:
             self.setup()
@@ -1091,12 +1070,12 @@ def add_param_from_top(self, param, target_component, setup=False):
         if target_component == "":
             setattr(self, param.name, param)
             self.top_level_params += [param.name]
-        else:
-            if target_component not in list(self.components.keys()):
-                raise AttributeError(
-                    "Can not find component '%s' in " "timing model." % target_component
-                )
+        elif target_component in list(self.components.keys()):
             self.components[target_component].add_param(param, setup=setup)
+        else:
+            raise AttributeError(
+                "Can not find component '%s' in " "timing model." % target_component
+            )
 
     def remove_param(self, param):
         """Remove a parameter from timing model.
@@ -1118,10 +1097,8 @@ def remove_param(self, param):
         self.setup()
 
     def get_params_mapping(self):
-        """Report whick component each parameter name comes from."""
-        param_mapping = {}
-        for p in self.top_level_params:
-            param_mapping[p] = "timing_model"
+        """Report which component each parameter name comes from."""
+        param_mapping = {p: "timing_model" for p in self.top_level_params}
         for cp in list(self.components.values()):
             for pp in cp.params:
                 param_mapping[pp] = cp.__class__.__name__
@@ -1218,10 +1195,9 @@ def delay(self, toas, cutoff_component="", include_last=True):
         Parameters
         ----------
         toas: pint.toa.TOAs
-            The toas for analysis delays.
+            The toas to analyze the delays.
         cutoff_component: str
-            The delay component name that a user wants the calculation to stop
-            at.
+            The delay component name at which to stop the calculation.
         include_last: bool
             If the cutoff delay component is included.
         """
@@ -1230,13 +1206,11 @@ def delay(self, toas, cutoff_component="", include_last=True):
             idx = len(self.DelayComponent_list)
         else:
             delay_names = [x.__class__.__name__ for x in self.DelayComponent_list]
-            if cutoff_component in delay_names:
-                idx = delay_names.index(cutoff_component)
-                if include_last:
-                    idx += 1
-            else:
+            if cutoff_component not in delay_names:
                 raise KeyError("No delay component named '%s'." % cutoff_component)
-
+            idx = delay_names.index(cutoff_component)
+            if include_last:
+                idx += 1
         # Do NOT cycle through delay_funcs - cycle through components until cutoff
         for dc in self.DelayComponent_list[:idx]:
             for df in dc.delay_funcs_component:
@@ -1259,30 +1233,26 @@ def phase(self, toas, abs_phase=None):
         # abs_phase defaults to True if AbsPhase is in the model, otherwise to
         # False.  Of course, if you manually set it, it will use that setting.
         if abs_phase is None:
-            if "AbsPhase" in list(self.components.keys()):
-                abs_phase = True
-            else:
-                abs_phase = False
+            abs_phase = "AbsPhase" in list(self.components.keys())
+        if not abs_phase:
+            return phase
         # If the absolute phase flag is on, use the TZR parameters to compute
         # the absolute phase.
-        if abs_phase:
-            if "AbsPhase" not in list(self.components.keys()):
-                # if no absolute phase (TZRMJD), add the component to the model and calculate it
-                from pint.models import absolute_phase
-
-                self.add_component(absolute_phase.AbsPhase(), validate=False)
-                self.make_TZR_toa(
-                    toas
-                )  # TODO:needs timfile to get all toas, but model doesn't have access to timfile. different place for this?
-                self.validate()
-            tz_toa = self.get_TZR_toa(toas)
-            tz_delay = self.delay(tz_toa)
-            tz_phase = Phase(np.zeros(len(toas.table)), np.zeros(len(toas.table)))
-            for pf in self.phase_funcs:
-                tz_phase += Phase(pf(tz_toa, tz_delay))
-            return phase - tz_phase
-        else:
-            return phase
+        if "AbsPhase" not in list(self.components.keys()):
+            # if no absolute phase (TZRMJD), add the component to the model and calculate it
+            from pint.models import absolute_phase
+
+            self.add_component(absolute_phase.AbsPhase(), validate=False)
+            self.make_TZR_toa(
+                toas
+            )  # TODO:needs timfile to get all toas, but model doesn't have access to timfile. different place for this?
+            self.validate()
+        tz_toa = self.get_TZR_toa(toas)
+        tz_delay = self.delay(tz_toa)
+        tz_phase = Phase(np.zeros(len(toas.table)), np.zeros(len(toas.table)))
+        for pf in self.phase_funcs:
+            tz_phase += Phase(pf(tz_toa, tz_delay))
+        return phase - tz_phase
 
     def total_dm(self, toas):
         """Calculate dispersion measure from all the dispersion type of components."""
@@ -1457,13 +1427,9 @@ def jump_flags_to_params(self, toas):
                 if tjv in tim_jump_values:
                     log.info(f"JUMP -tim_jump {tjv} already exists")
                     tim_jump_values.remove(tjv)
-        if used_indices:
-            num = max(used_indices) + 1
-        else:
-            num = 1
-
+        num = max(used_indices) + 1 if used_indices else 1
         if not tim_jump_values:
-            log.info(f"All tim_jump values have corresponding JUMPs")
+            log.info("All tim_jump values have corresponding JUMPs")
             return
 
         # FIXME: arrange for these to be in a sensible order (might not be integers
@@ -1509,7 +1475,7 @@ def delete_jump_and_flags(self, toa_table, jump_num):
             Specifies the index of the jump to be deleted.
         """
         # remove jump of specified index
-        self.remove_param("JUMP" + str(jump_num))
+        self.remove_param(f"JUMP{str(jump_num)}")
 
         # remove jump flags from selected TOA tables
         if toa_table is not None:
@@ -1736,10 +1702,7 @@ def d_phase_d_param_num(self, toas, param, step=1e-2):
         par = getattr(self, param)
         ori_value = par.value
         unit = par.units
-        if ori_value == 0:
-            h = 1.0 * step
-        else:
-            h = ori_value * step
+        h = 1.0 * step if ori_value == 0 else ori_value * step
         parv = [par.value - h, par.value + h]
 
         phase_i = (
@@ -1773,10 +1736,7 @@ def d_delay_d_param_num(self, toas, param, step=1e-2):
             log.warning("Parameter '%s' is not used by timing model." % param)
             return np.zeros(toas.ntoas) * (u.second / par.units)
         unit = par.units
-        if ori_value == 0:
-            h = 1.0 * step
-        else:
-            h = ori_value * step
+        h = 1.0 * step if ori_value == 0 else ori_value * step
         parv = [par.value - h, par.value + h]
         delay = np.zeros((toas.ntoas, 2))
         for ii, val in enumerate(parv):
diff --git a/src/pint/simulation.py b/src/pint/simulation.py
index 528b77865..d3602504d 100644
--- a/src/pint/simulation.py
+++ b/src/pint/simulation.py
@@ -11,7 +11,7 @@
 import pint.toa
 import pint.models
 from pint.residuals import Residuals
-from pint.observatory import Observatory, bipm_default, get_observatory
+from pint.observatory import bipm_default, get_observatory
 
 __all__ = [
     "zero_residuals",
@@ -65,10 +65,7 @@ def zero_residuals(ts, model, maxiter=10, tolerance=None):
     ts.compute_pulse_numbers(model)
     maxresid = None
     if tolerance is None:
-        if pint.utils.check_longdouble_precision():
-            tolerance = 1 * u.ns
-        else:
-            tolerance = 5 * u.us
+        tolerance = 1 * u.ns if pint.utils.check_longdouble_precision() else 5 * u.us
     for i in range(maxiter):
         r = Residuals(ts, model, track_mode="use_pulse_numbers")
         resids = r.calc_time_resids(calctype="taylor")