feature: solar wind proxy regression

CHANGELOG-unreleased.md

@@ -14,6 +14,8 @@ the released changes.
 - `pintk` Diff/Unc calculation now uses post-fit uncertainties.
 ### Added
 - Plot whitened DM residuals in pintk.
+- Time-domain solar wind GP noise components: ridge, squared-exponential, Matérn, and quasi-periodic kernels
+- Regression tests for noise design-matrix caching, including multi-basis coverage (red, DMGP, SWGP, and chromatic GP)
 ### Fixed
 - `WidebandTOAFitter` raises a warning if the model has correlated errors (It used to give wrong results before).
 - Fixed bug where "include_bipm" flag was being ignored when loading Fermi TOAs with weights, now defaults to using EPHEM, CLOCK and PLANET_SHAPIRO from the timing model

src/pint/fitter.py

@@ -94,7 +94,7 @@
 )
 from pint.residuals import Residuals, WidebandTOAResiduals
 from pint.toa import TOAs
-from pint.utils import FTest, normalize_designmatrix
+from pint.utils import FTest, get_phiinv, normalize_designmatrix
 
 __all__ = [
     "Fitter",
@@ -1341,7 +1341,9 @@ def step(self) -> np.ndarray:
         else:
             M, params, units = self.model.full_designmatrix(self.fitter.toas)
             M, norm = normalize_designmatrix(M, params)
-            phiinv = 1 / self.model.full_basis_weight(self.fitter.toas) / norm**2
+            phi = self.model.full_basis_weight(self.fitter.toas)
+            phiinv = get_phiinv(phi)
+            phiinv = (phiinv / norm).T / norm  # normalize the phi^-1 matrix
             Nvec = (
                 self.model.scaled_toa_uncertainty(self.fitter.toas).to_value(u.s) ** 2
             )
@@ -1902,7 +1904,9 @@ def fit_toas(
             else:
                 M, params, units = self.model.full_designmatrix(self.toas)
                 M, norm = normalize_designmatrix(M, params)
-                phiinv = 1 / self.model.full_basis_weight(self.toas) / norm**2
+                phi = self.model.full_basis_weight(self.toas)
+                phiinv = get_phiinv(phi)
+                phiinv = (phiinv / norm).T / norm  # normalize the phi^-1 matrix
                 Nvec = self.model.scaled_toa_uncertainty(self.toas).to(u.s).value ** 2
                 mtcm, mtcy = get_gls_mtcm_mtcy(phiinv, Nvec, M, residuals)
 
@@ -2205,13 +2209,19 @@ def fit_toas(
                 phi = self.model.noise_model_basis_weight(self.toas)
                 phiinv = np.zeros(M.shape[1])
                 if Mn is not None and phi is not None:
-                    phiinv = np.concatenate((phiinv, 1 / phi))
+                    n_tm = M.shape[1]
                     new_d_matrix = combine_design_matrices_by_param(d_matrix, Mn)
                     M, params, units = (
                         new_d_matrix.matrix,
                         new_d_matrix.derivative_params,
                         new_d_matrix.param_units,
                     )
+                    if np.ndim(phi) == 1:
+                        phiinv = np.zeros(M.shape[1])
+                        phiinv[n_tm:] = 1 / phi
+                    else:
+                        phiinv = np.zeros((M.shape[1], M.shape[1]))
+                        phiinv[n_tm:, n_tm:] = get_phiinv(phi)
 
             ntmpar = self.model.ntmpar
 
@@ -2228,11 +2238,17 @@ def fit_toas(
                 mtcy = np.dot(cm.T, residuals)
                 # mtcm, mtcy = get_gls_mtcm_mtcy_fullcov(cov, M, residuals)
             else:
-                phiinv /= norm**2
+                if np.ndim(phiinv) == 1:
+                    phiinv /= norm**2
+                else:
+                    phiinv = (phiinv / norm).T / norm
                 Nvec = self.scaled_all_sigma() ** 2
                 cinv = 1 / Nvec
                 mtcm = np.dot(M.T, cinv[:, None] * M)
-                mtcm += np.diag(phiinv)
+                if np.ndim(phiinv) == 1:
+                    mtcm += np.diag(phiinv)
+                else:
+                    mtcm += phiinv
                 mtcy = np.dot(M.T, cinv * residuals)
                 # mtcm, mtcy = get_gls_mtcm_mtcy(phiinv, Nvec, M, residuals)
 
@@ -2249,7 +2265,11 @@ def fit_toas(
             if full_cov:
                 chi2 = np.dot(newres, scipy.linalg.cho_solve(cf, newres))
             else:
-                chi2 = np.dot(newres, cinv * newres) + np.dot(xhat, phiinv * xhat)
+                if np.ndim(phiinv) == 1:
+                    prior_chi2 = np.dot(xhat, phiinv * xhat)
+                else:
+                    prior_chi2 = np.dot(xhat, phiinv @ xhat)
+                chi2 = np.dot(newres, cinv * newres) + prior_chi2
 
             # compute absolute estimates, normalized errors, covariance matrix
             dpars = xhat / norm
@@ -2631,7 +2651,10 @@ def get_gls_mtcm_mtcy(
     """
     cinv = 1 / Nvec
     mtcmplain = np.dot(M.T, cinv[:, None] * M)
-    mtcm = mtcmplain + np.diag(phiinv)
+    if np.ndim(phiinv) == 1:
+        mtcm = mtcmplain + np.diag(phiinv)
+    else:
+        mtcm = mtcmplain + phiinv
     mtcy = np.dot(M.T, cinv * residuals)
     return (mtcm, mtcy) if not return_mtcmplain else (mtcm, mtcy, mtcmplain)
 

src/pint/models/__init__.py

@@ -49,6 +49,7 @@
     PLRedNoise,
     ScaleToaError,
     ScaleDmError,
+    TimeDomainSWNoise,
 )
 from pint.models.phase_offset import PhaseOffset
 from pint.models.piecewise import PiecewiseSpindown

src/pint/models/model_builder.py

@@ -274,6 +274,13 @@ def _validate_components(self):
                 if v.category == "pulsar_system":
                     # The pulsar system will be selected by parameter BINARY
                     continue
+                elif self.all_components.components[superset].category == v.category:
+                    # Same-category components are mutually exclusive alternatives.
+                    # A superset component in the same category is selected over a
+                    # subset when its extra parameters appear in the parfile, so
+                    # there is no ambiguity.  (E.g. SolarWindProxyRegression vs
+                    # SolarWindDispersion, discriminated by SWPRBETA1.)
+                    continue
                 else:
                     raise ComponentConflict(m)