search

nanograv / PINT

PINT is not TEMPO3 -- Software for high-precision pulsar timing
Other
118 stars 101 forks source link

model.toa_covariance_matrix is wrong if noise components are present #978

Closed aarchiba closed 3 years ago

aarchiba commented 3 years ago

In PINT's toa_covariance_matrix, if there is no noise model, you get a diagonal matrix with the square of each TOA's uncertainty on the diagonal. If there is a noise model, though, this term is not added, only the covariances coming from the noise components - so in particular if your noise model is just ECORRs the covariance matrix has zero eigenvalues and Cholesky decomposition fails. I do not think this can be correct.

https://github.com/nanograv/PINT/blob/9250b5a1d487acf866bdd5605ebd212d758bbb4e/src/pint/models/timing_model.py#L1106

import io
import re
from copy import deepcopy

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
import pytest
from astropy import log
from astropy.time import TimeDelta
from scipy.linalg import block_diag, cho_factor, cho_solve, cholesky, eigvalsh

import pint.fitter
from pint.models import get_model
from pint.models.timing_model import MissingTOAs
from pint.toa import make_fake_toas, merge_TOAs

par_eccentric = """
PSR J1234+5678
F0 500 0 1
ELAT 0 0 1
ELONG 0 0 1
PEPOCH 57000
POSEPOCH 57000
DM 10 0 1
SOLARN0 0 0 1
BINARY BT
PB 1 0 1
A1 10 0 1
ECC 0.95 0 1
OM 0 0 1
T0 57000 0 1
"""
model_eccentric = get_model(io.StringIO(par_eccentric))

toas_ecorr = merge_TOAs(
    [
        make_fake_toas(57000, 57001, 100, model_eccentric, freq=1400 * u.MHz, obs="@"),
        make_fake_toas(57000, 57001, 100, model_eccentric, freq=1800 * u.MHz, obs="@"),
    ]
)
toas_ecorr.adjust_TOAs(TimeDelta(np.random.randn(len(toas_ecorr)) * toas_ecorr.table["error"]))

model_wrong_ecorr = get_model(io.StringIO("\n".join([par_eccentric, "ECORR tel @ 1"])))
model_wrong_ecorr.ECC.value = 0.5

f6 = pint.fitter.GLSFitter(toas_ecorr, model_wrong_ecorr)
f6.model.free_params = ["ECC"]
f6.fit_toas(full_cov=True)

yields:

DEBUG: Skipping Shapiro delay for Barycentric TOAs [pint.models.solar_system_shapiro]
DEBUG: Skipping Shapiro delay for Barycentric TOAs [pint.models.solar_system_shapiro]
DEBUG: Skipping Shapiro delay for Barycentric TOAs [pint.models.solar_system_shapiro]
DEBUG: Skipping Shapiro delay for Barycentric TOAs [pint.models.solar_system_shapiro]
DEBUG: Skipping Shapiro delay for Barycentric TOAs [pint.models.solar_system_shapiro]
DEBUG: Skipping Shapiro delay for Barycentric TOAs [pint.models.solar_system_shapiro]
DEBUG: Skipping Shapiro delay for Barycentric TOAs [pint.models.solar_system_shapiro]
DEBUG: No pulse numbers found in the TOAs [pint.toa]
INFO: Applying clock corrections (include_gps = True, include_bipm = False) [pint.toa]
INFO: Special observatory location. No clock corrections applied. [pint.observatory.special_locations]
INFO: Computing TDB columns. [pint.toa]
INFO: Using EPHEM = DE421 for TDB calculation. [pint.toa]
INFO: Computing PosVels of observatories and Earth, using DE421 [pint.toa]
DEBUG: SSB obs pos [0. 0. 0.] m [pint.toa]
DEBUG: Set solar system ephemeris to (cached) link https://data.nanograv.org/static/data/ephem/de421.bsp [pint.solar_system_ephemerides]
DEBUG: Adding columns ssb_obs_pos ssb_obs_vel obs_sun_pos [pint.toa]
DEBUG: Skipping Shapiro delay for Barycentric TOAs [pint.models.solar_system_shapiro]
DEBUG: Adding column ssb_obs_vel_ecl [pint.toa]
/home/archibald/.virtualenvs/pint-dev3/lib/python3.8/site-packages/astropy/table/row.py:76: FutureWarning: elementwise == comparison failed and returning scalar instead; this will raise an error or perform elementwise comparison in the future.
  return self.as_void() == other
---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
<ipython-input-81-e9472e414c6c> in <module>
      1 f6 = pint.fitter.GLSFitter(toas_ecorr, model_wrong_ecorr)
      2 f6.model.free_params = ["ECC"]
----> 3 f6.fit_toas(full_cov=True)

~/projects/pint/PINT/src/pint/fitter.py in fit_toas(self, maxiter, threshold, full_cov)
   1539             if full_cov:
   1540                 cov = self.model.toa_covariance_matrix(self.toas)
-> 1541                 cf = sl.cho_factor(cov)
   1542                 cm = sl.cho_solve(cf, M)
   1543                 mtcm = np.dot(M.T, cm)

~/.virtualenvs/pint-dev3/lib/python3.8/site-packages/scipy/linalg/decomp_cholesky.py in cho_factor(a, lower, overwrite_a, check_finite)
    150 
    151     """
--> 152     c, lower = _cholesky(a, lower=lower, overwrite_a=overwrite_a, clean=False,
    153                          check_finite=check_finite)
    154     return c, lower

~/.virtualenvs/pint-dev3/lib/python3.8/site-packages/scipy/linalg/decomp_cholesky.py in _cholesky(a, lower, overwrite_a, clean, check_finite)
     35     c, info = potrf(a1, lower=lower, overwrite_a=overwrite_a, clean=clean)
     36     if info > 0:
---> 37         raise LinAlgError("%d-th leading minor of the array is not positive "
     38                           "definite" % info)
     39     if info < 0:

LinAlgError: 101-th leading minor of the array is not positive definite

Also:

eigvalsh(model_wrong_ecorr.toa_covariance_matrix(toas_ecorr))

yields

array([0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00,
       0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00,
       0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00,
       0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00,
       0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00,
       0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00,
       0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00,
       0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00,
       0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00,
       0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00,
       0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00,
       0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00, 0.e+00,
       0.e+00, 0.e+00, 0.e+00, 0.e+00, 2.e-12, 2.e-12, 2.e-12, 2.e-12,
       2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12,
       2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12,
       2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12,
       2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12,
       2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12,
       2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12,
       2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12,
       2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12,
       2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12,
       2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12,
       2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12, 2.e-12,
aarchiba commented 3 years ago

Fixing it seems to lead to discrepancies with tempo2 and with the reduced-rank versions, although not large discrepancies. Test failures, though.

aarchiba commented 3 years ago

The problem occurs specifically if you have a noise model but not EFAC/EQUAD.