Fixed bug in residual calculation when adding or removing phase wraps…

… when not fully phase-connected. Also cleaned up residual calculation.

src/pint/pintk/pulsar.py

@@ -5,6 +5,7 @@
 pre/post fit residuals, and other useful information.
 self.selected_toas = selected toas, self.all_toas = all toas in tim file
 """
+
 import copy
 
 import astropy.units as u
@@ -100,8 +101,9 @@ def __init__(self, parfile=None, timfile=None, ephem=None, fitter="GLSFitter"):
 
         self.all_toas.print_summary()
 
-        self.prefit_resids = Residuals(self.all_toas, self.prefit_model)
-        self.selected_prefit_resids = self.prefit_resids
+        self.use_pulse_numbers = False
+        self.fitted = False
+        self.update_resids()
         print(
             "RMS pre-fit PINT residuals are %.3f us\n"
             % self.prefit_resids.rms_weighted().to(u.us).value
@@ -121,11 +123,9 @@ def __init__(self, parfile=None, timfile=None, ephem=None, fitter="GLSFitter"):
         else:
             self.fit_method = fitter
         self.fitter = None
-        self.fitted = False
         self.stashed = None  # for temporarily stashing some TOAs
         self.faketoas1 = None  # for random models
         self.faketoas = None  # for random models
-        self.use_pulse_numbers = False
 
     @property
     def name(self):
@@ -205,18 +205,36 @@ def update_resids(self):
         # update the pre and post fit residuals using all_toas
         track_mode = "use_pulse_numbers" if self.use_pulse_numbers else None
         self.prefit_resids = Residuals(
-            self.all_toas, self.prefit_model, track_mode=track_mode
+            self.all_toas, self.prefit_model, subtract_mean=False, track_mode=track_mode
         )
         if self.selected_toas.ntoas and self.selected_toas.ntoas != self.all_toas.ntoas:
             self.selected_prefit_resids = Residuals(
-                self.selected_toas, self.prefit_model, track_mode=track_mode
+                self.selected_toas,
+                self.prefit_model,
+                subtract_mean=False,
+                track_mode=track_mode,
             )
         else:
             self.selected_prefit_resids = self.prefit_resids
         if self.fitted:
             self.postfit_resids = Residuals(
-                self.all_toas, self.postfit_model, track_mode=track_mode
+                self.all_toas,
+                self.postfit_model,
+                subtract_mean=False,
+                track_mode=track_mode,
             )
+            if (
+                self.selected_toas.ntoas
+                and self.selected_toas.ntoas != self.all_toas.ntoas
+            ):
+                self.selected_postfit_resids = Residuals(
+                    self.selected_toas,
+                    self.postfit_model,
+                    subtract_mean=False,
+                    track_mode=track_mode,
+                )
+            else:
+                self.selected_postfit_resids = self.postfit_resids
 
     def orbitalphase(self):
         """
@@ -470,7 +488,7 @@ def print_chi2(self, selected):
             self.prefit_resids = self.postfit_resids
             self.add_model_params()
 
-        self.selected_resids = Residuals(self.selected_toas, self.prefit_model)
+        self.update_resids()
 
         wrms = self.selected_resids.rms_weighted()
         print("------------------------------------")
@@ -496,12 +514,7 @@ def fit(self, selected, iters=4, compute_random=False):
             self.prefit_resids = self.postfit_resids
             self.add_model_params()
 
-        if self.selected_toas.ntoas != self.all_toas.ntoas:
-            self.selected_prefit_resids = Residuals(
-                self.selected_toas, self.prefit_model
-            )
-        else:
-            self.selected_prefit_resids = self.prefit_resids
+        self.update_resids()
 
         # Have to change the fitter for each fit since TOAs and models change
         log.info(f"Using {self.fit_method}")
@@ -533,12 +546,7 @@ def fit(self, selected, iters=4, compute_random=False):
         self.selected_toas.compute_pulse_numbers(self.postfit_model)
 
         # Compute the residuals using correct pulse numbers
-        self.postfit_resids = Residuals(self.all_toas, self.postfit_model)
-        self.selected_postfit_resids = (
-            self.postfit_resids
-            if np.all(selected)
-            else Residuals(self.selected_toas, self.postfit_model)
-        )
+        self.update_resids()
 
         # Need this since it isn't updated using self.fitter.update_model()
         self.fitter.model.CHI2.value = self.selected_postfit_resids.chi2
@@ -572,6 +580,8 @@ def fit(self, selected, iters=4, compute_random=False):
                 getattr(pm_no_jumps, param).value = 0.0
                 getattr(pm_no_jumps, param).frozen = True
         self.prefit_resids_no_jumps = Residuals(self.all_toas, pm_no_jumps)
+        self.update_resids()
+        self.prefit_resids_no_jumps = self.prefit_resids
 
         # Store some key params for possible F-testing
         self.lastfit = {