diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 8200caa75..32b491913 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -29,6 +29,14 @@ the released changes. - Piecewise orbital model (`BinaryBTPiecewise`) - `TimingModel.fittable_params` property - Simulate correlated noise using `pint.simulation` (also available via the `zima` script) +- Optionally return the the log normalization factor of the likelihood function from the `Residuals.calc_chi2()` method. +- `DownhilWLSFitter` can now estimate white noise parameters and their uncertainties. +- `Residuals.lnlikelihood()` method +- `pint.utils.akaike_information_criterion()` function +- `TimingModel.d_toasigma_d_param` method to compute derivatives of scaled TOA uncertainties w.r.t. white noise parameters. +- `TimingModel.toasigma_derivs` property to get all derivatives functions of scaled TOA uncertainties. +- `ScaleToaError.register_toasigma_deriv_funcs` method to populate derivatives of scaled TOA uncertainties. +- `ScaleToaError.d_toasigma_d_EFAC` and `ScaleToaError.d_toasigma_d_EQUAD` methods. ### Fixed - Wave model `validate()` can correctly use PEPOCH to assign WAVEEPOCH parameter - Fixed RTD by specifying theme explicitly. @@ -39,4 +47,5 @@ the released changes. - `get_TOAs` raises an exception upon finding mixed narrowband and wideband TOAs in a tim file. `TOAs.is_wideband` returns True only if *ALL* TOAs have the -pp_dm flag. - `TimingModel.designmatrix()` method will fail with an informative error message if there are free unfittable parameters in the timing model. - `make_fake_toas_uniform` and `make_fake_toas_fromMJDs` respects units of errors +- Robust access of EPHEM and PLANET_SHAPIRO in `make_fake_toas_fromtim` ### Removed diff --git a/requirements.txt b/requirements.txt index 1c483a0a8..16f76401e 100644 --- a/requirements.txt +++ b/requirements.txt @@ -9,3 +9,4 @@ corner>=2.0.1 uncertainties loguru nestle>=0.2.0 +numdifftools diff --git a/src/pint/fitter.py b/src/pint/fitter.py index ce678edc2..287f7f467 100644 --- a/src/pint/fitter.py +++ b/src/pint/fitter.py @@ -67,6 +67,7 @@ import scipy.linalg import scipy.optimize as opt from loguru import logger as log +from numdifftools import Hessdiag import pint import pint.utils @@ -1166,7 +1167,7 @@ def __init__(self, toas, model, track_mode=None, residuals=None): ) self.method = "downhill_checked" - def fit_toas( + def _fit_toas( self, maxiter=20, required_chi2_decrease=1e-2, @@ -1174,37 +1175,11 @@ def fit_toas( min_lambda=1e-3, debug=False, ): - """Carry out a cautious downhill fit. - - This tries to take the same steps as - :func:`pint.fitter.WLSFitter.fit_toas` or - :func:`pint.fitter.GLSFitter.fit_toas` or - :func:`pint.fitter.WidebandTOAFitter.fit_toas`. At each step, it - checks whether the new model has a better ``chi2`` than the current - one; if the new model is invalid or worse than the current one, it - tries taking a shorter step in the same direction. This can exit if it - exceeds the maximum number of iterations or if improvement is not - possible even with very short steps, or it can exit successfully if a - full-size step is taken and it does not decrease the ``chi2`` by much. + """Downhill fit implementation for fitting the timing model parameters. + The `fit_toas()` calls this method iteratively to fit the timing model parameters + while also fitting for white noise parameters. - The attribute ``self.converged`` is set to True or False depending on - whether the process actually converged. - - Parameters - ========== - - maxiter : int - Abandon the process if this many successful steps have been taken. - required_chi2_decrease : float - A full-size step that makes less than this much improvement is taken - to indicate that the fitter has converged. - max_chi2_increase : float - If this is positive, consider taking steps that slightly worsen the chi2 in hopes - of eventually finding our way downhill. - min_lambda : float - If steps are shrunk by this factor and still don't result in improvement, abandon hope - of convergence and stop. - """ + See documentation of the `fit_toas()` method for more details.""" # setup self.model.validate() self.model.validate_toas(self.toas) @@ -1213,6 +1188,7 @@ def fit_toas( self.converged = False # algorithm exception = None + for i in range(maxiter): step = current_state.step lambda_ = 1 @@ -1266,6 +1242,7 @@ def fit_toas( log.debug( f"Stopping because maximum number of iterations ({maxiter}) reached" ) + self.current_state = best_state # collect results self.model = self.current_state.model @@ -1277,6 +1254,7 @@ def fit_toas( self.parameter_correlation_matrix = ( self.parameter_covariance_matrix.to_correlation_matrix() ) + for p, e in zip(self.current_state.params, self.errors): try: # I don't know why this fails with multiprocessing, but bypass if it does @@ -1297,10 +1275,154 @@ def fit_toas( raise MaxiterReached(f"Convergence not detected after {maxiter} steps.") return self.converged + def fit_toas( + self, + maxiter=20, + noise_fit_niter=5, + required_chi2_decrease=1e-2, + max_chi2_increase=1e-2, + min_lambda=1e-3, + noisefit_method="Newton-CG", + debug=False, + ): + """Carry out a cautious downhill fit. + + This tries to take the same steps as + :func:`pint.fitter.WLSFitter.fit_toas` or + :func:`pint.fitter.GLSFitter.fit_toas` or + :func:`pint.fitter.WidebandTOAFitter.fit_toas`. At each step, it + checks whether the new model has a better ``chi2`` than the current + one; if the new model is invalid or worse than the current one, it + tries taking a shorter step in the same direction. This can exit if it + exceeds the maximum number of iterations or if improvement is not + possible even with very short steps, or it can exit successfully if a + full-size step is taken and it does not decrease the ``chi2`` by much. + + The attribute ``self.converged`` is set to True or False depending on + whether the process actually converged. + + This function can also estimate white noise parameters (EFACs and EQUADs) + and their uncertainties. + + If there are no free white noise parameters, this function will do one + iteration of the downhill fit (implemented in the `_fit_toas()` method). + If free white noise parameters are present, it will fit for them by numerically + maximizing the likelihood function (implemented in the `_fit_noise()` method). + The timing model fit and the noise model fit are run iteratively in an alternating + fashion. Fitting for a white noise parameter is as simple as:: + + fitter.model.EFAC1.frozen = False + fitter.fit_toas() + + + Parameters + ========== + + maxiter : int + Abandon the process if this many successful steps have been taken. + required_chi2_decrease : float + A full-size step that makes less than this much improvement is taken + to indicate that the fitter has converged. + max_chi2_increase : float + If this is positive, consider taking steps that slightly worsen the chi2 in hopes + of eventually finding our way downhill. + min_lambda : float + If steps are shrunk by this factor and still don't result in improvement, abandon hope + of convergence and stop. + noisefit_method: str + Algorithm used to fit for noise parameters. See the documentation for + `scipy.optimize.minimize()` for more details and available options. + """ + free_noise_params = self._get_free_noise_params() + + if len(free_noise_params) == 0: + return self._fit_toas( + maxiter=maxiter, + required_chi2_decrease=required_chi2_decrease, + max_chi2_increase=required_chi2_decrease, + min_lambda=required_chi2_decrease, + debug=debug, + ) + else: + log.debug("Will fit for noise parameters.") + for _ in range(noise_fit_niter): + self._fit_toas( + maxiter=maxiter, + required_chi2_decrease=required_chi2_decrease, + max_chi2_increase=max_chi2_increase, + min_lambda=min_lambda, + debug=debug, + ) + values, errors = self._fit_noise(noisefit_method=noisefit_method) + self._update_noise_params(values, errors) + + return self._fit_toas( + maxiter=maxiter, + required_chi2_decrease=required_chi2_decrease, + max_chi2_increase=max_chi2_increase, + min_lambda=min_lambda, + debug=debug, + ) + @property def fac(self): return self.current_state.fac + def _get_free_noise_params(self): + """Returns a list of all free noise parameters.""" + return [ + fp + for fp in self.model.get_params_of_component_type("NoiseComponent") + if not getattr(self.model, fp).frozen + ] + + def _update_noise_params(self, values, errors): + """Update the model using estimated noise parameters.""" + free_noise_params = self._get_free_noise_params() + for fp, val, err in zip(free_noise_params, values, errors): + getattr(self.model, fp).value = val + getattr(self.model, fp).uncertainty_value = err + + def _fit_noise(self, noisefit_method="Newton-CG"): + """Estimate noise parameters and their uncertainties. Noise parameters + are estimated by numerically maximizing the log-likelihood function including + the normalization term. The uncertainties thereof are computed using the + numerically-evaluated Hessian.""" + free_noise_params = self._get_free_noise_params() + + xs0 = [getattr(self.model, fp).value for fp in free_noise_params] + + model1 = copy.deepcopy(self.model) + res = Residuals(self.toas, model1) + + def _mloglike(xs): + """Negative of the log-likelihood function.""" + for fp, x in zip(free_noise_params, xs): + getattr(res.model, fp).value = x + + return -res.lnlikelihood() + + def _mloglike_grad(xs): + """Gradient of the negative of the log-likelihood function w.r.t. white noise parameters.""" + for fp, x in zip(free_noise_params, xs): + getattr(res.model, fp).value = x + + return np.array( + [ + -res.d_lnlikelihood_d_whitenoise_param(par).value + for par in free_noise_params + ] + ) + + maxlike_result = opt.minimize( + _mloglike, xs0, method=noisefit_method, jac=_mloglike_grad + ) + + hess = Hessdiag(_mloglike) + errs = np.sqrt(1 / hess(maxlike_result.x)) + + return maxlike_result.x, errs + class WLSState(ModelState): def __init__(self, fitter, model, threshold=None): diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index 49036fd8a..095d6ffab 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -58,7 +58,7 @@ def __init__( name="EFAC", units="", aliases=["T2EFAC", "TNEF"], - description="A multiplication factor for the measured TOA uncertainties,", + description="A multiplication factor on the measured TOA uncertainties,", ) ) @@ -67,7 +67,7 @@ def __init__( name="EQUAD", units="us", aliases=["T2EQUAD"], - description="An error term added in quadrature to the TOA uncertainty.", + description="An error term added in quadrature to the scaled (by EFAC) TOA uncertainty.", ) ) @@ -75,11 +75,12 @@ def __init__( maskParameter( name="TNEQ", units=u.LogUnit(physical_unit=u.second), - description="A log10-scale error term added in quadrature to the TOA uncertainty", + description="An error term added in quadrature to the scaled (by EFAC) TOA uncertainty in units of log10(second).", ) ) self.covariance_matrix_funcs += [self.sigma_scaled_cov_matrix] self.scaled_toa_sigma_funcs += [self.scale_toa_sigma] + self.toasigma_deriv_funcs = {} def setup(self): super().setup() @@ -117,9 +118,7 @@ def setup(self): units="us", index=tneq_par.index, aliases=["T2EQUAD"], - description="An error term " - " added in quadrature to the" - " scaled (by EFAC) TOA uncertainty.", + description="An error term added in quadrature to the scaled (by EFAC) TOA uncertainty.", ) ) EQUAD_par = getattr(self, EQUAD_name) @@ -131,6 +130,21 @@ def setup(self): par = getattr(self, pp) self.EQUADs[pp] = (par.key, par.key_value) + for ef in self.EFACs: + self.register_toasigma_deriv_funcs(self.d_toasigma_d_EFAC, ef) + + for eq in self.EQUADs: + self.register_toasigma_deriv_funcs(self.d_toasigma_d_EQUAD, eq) + + def register_toasigma_deriv_funcs(self, func, param): + pn = self.match_param_aliases(param) + if pn not in list(self.toasigma_deriv_funcs.keys()): + self.toasigma_deriv_funcs[pn] = [func] + elif func in self.toasigma_deriv_funcs[pn]: + return + else: + self.toasigma_deriv_funcs[pn] += [func] + def validate(self): super().validate() # check duplicate @@ -139,7 +153,7 @@ def validate(self): if [x for x in l if l.count(x) > 1] != []: raise ValueError(f"'{el}' have duplicated keys and key values.") - def scale_toa_sigma(self, toas): + def scale_toa_sigma(self, toas, warn=True): sigma_scaled = toas.table["error"].quantity.copy() for equad_name in self.EQUADs: equad = getattr(self, equad_name) @@ -148,14 +162,14 @@ def scale_toa_sigma(self, toas): mask = equad.select_toa_mask(toas) if np.any(mask): sigma_scaled[mask] = np.hypot(sigma_scaled[mask], equad.quantity) - else: + elif warn: warnings.warn(f"EQUAD {equad} has no TOAs") for efac_name in self.EFACs: efac = getattr(self, efac_name) mask = efac.select_toa_mask(toas) if np.any(mask): sigma_scaled[mask] *= efac.quantity - else: + elif warn: warnings.warn(f"EFAC {efac} has no TOAs") return sigma_scaled @@ -163,6 +177,39 @@ def sigma_scaled_cov_matrix(self, toas): scaled_sigma = self.scale_toa_sigma(toas).to(u.s).value ** 2 return np.diag(scaled_sigma) + def d_toasigma_d_EFAC(self, toas, param): + par = getattr(self, param) + mask = par.select_toa_mask(toas) + result = np.zeros(len(toas)) << u.s + result[mask] = self.scale_toa_sigma(toas[mask], warn=False).to( + u.s + ) / par.quantity.to(u.dimensionless_unscaled) + return result + + def d_toasigma_d_EQUAD(self, toas, param): + par = getattr(self, param) + mask = par.select_toa_mask(toas) + toas_mask = toas[mask] + + result = np.zeros(len(toas)) << u.dimensionless_unscaled + + sigma_mask = self.scale_toa_sigma(toas_mask, warn=False) + + sigma2_mask_noefac = toas_mask.get_errors().to(u.s) ** 2 + for equad_name in self.EQUADs: + equad = getattr(self, equad_name) + if equad.quantity is None: + continue + eqmask = equad.select_toa_mask(toas_mask) + if np.any(eqmask): + sigma2_mask_noefac[eqmask] += equad.quantity**2 + + result[mask] = (sigma_mask * par.quantity / sigma2_mask_noefac).to( + u.dimensionless_unscaled + ) + + return result + class ScaleDmError(NoiseComponent): """Correction for estimated wideband DM measurement uncertainty. @@ -190,8 +237,7 @@ def __init__( maskParameter( name="DMEFAC", units="", - description="A multiplication factor on" - " the measured DM uncertainties,", + description="A multiplication factor on the measured DM uncertainties,", ) ) @@ -199,9 +245,7 @@ def __init__( maskParameter( name="DMEQUAD", units="pc / cm ^ 3", - description="An error term added in " - "quadrature to the scaled (by" - " EFAC) TOA uncertainty.", + description="An error term added in quadrature to the scaled (by EFAC) TOA uncertainty.", ) ) @@ -428,7 +472,7 @@ def __init__( name="TNDMAMP", units="", aliases=[], - description="Amplitude of powerlaw " "DM noise in tempo2 format", + description="Amplitude of powerlaw DM noise in tempo2 format", ) ) self.add_param( @@ -537,7 +581,7 @@ def __init__( name="RNAMP", units="", aliases=[], - description="Amplitude of powerlaw " "red noise.", + description="Amplitude of powerlaw red noise.", ) ) self.add_param( @@ -545,7 +589,7 @@ def __init__( name="RNIDX", units="", aliases=[], - description="Spectral index of " "powerlaw red noise.", + description="Spectral index of powerlaw red noise.", ) ) @@ -554,7 +598,7 @@ def __init__( name="TNREDAMP", units="", aliases=[], - description="Amplitude of powerlaw " "red noise in tempo2 format", + description="Amplitude of powerlaw red noise in tempo2 format", ) ) self.add_param( @@ -562,7 +606,7 @@ def __init__( name="TNREDGAM", units="", aliases=[], - description="Spectral index of powerlaw " "red noise in tempo2 format", + description="Spectral index of powerlaw red noise in tempo2 format", ) ) self.add_param( diff --git a/src/pint/models/parameter.py b/src/pint/models/parameter.py index 03daac7a6..c94e62f6e 100644 --- a/src/pint/models/parameter.py +++ b/src/pint/models/parameter.py @@ -158,6 +158,8 @@ class Parameter: use_alias : str or None Alias to use on write; normally whatever alias was in the par file it was read from + parent: pint.models.timing_model.Component, optional + The parent timing model component Attributes ---------- @@ -183,6 +185,7 @@ def __init__( # The input parameter from parfile, which can be an alias of the parameter # TODO give a better name and make it easy to access. self._parfile_name = name + self.units = units # Default unit self.quantity = value # The value of parameter, internal storage self.prior = prior @@ -216,6 +219,7 @@ def quantity(self, val): raise ValueError("Setting an existing value to None is not allowed.") self._quantity = val return + self._quantity = self._set_quantity(val) @property @@ -240,6 +244,7 @@ def value(self, val): ) else: self.value = val + self._quantity = self._set_quantity(val) @property @@ -969,8 +974,8 @@ class intParameter(Parameter): ---------- name : str The name of the parameter. - value : str, bool, [0,1] - The input parameter boolean value. + value : int + The parameter value. description : str, optional A short description of what this parameter means. aliases : list, optional diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 0b6350548..20eb85aba 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -625,7 +625,16 @@ def fittable_params(self): return [ p for p in self.params - if (p in self.phase_deriv_funcs or p in self.delay_deriv_funcs) + if ( + p in self.phase_deriv_funcs + or p in self.delay_deriv_funcs + or ( + ( + hasattr(self, "toasigma_deriv_funcs") + and p in self.toasigma_deriv_funcs + ) + ) + ) ] def match_param_aliases(self, alias): @@ -1058,9 +1067,14 @@ def delay_deriv_funcs(self): @property_exists def dm_derivs(self): # TODO need to be careful about the name here. - """List of dm derivative functions.""" + """List of DM derivative functions.""" return self.get_deriv_funcs("DelayComponent", "dm") + @property_exists + def toasigma_derivs(self): + """List of scaled TOA uncertainty derivative functions""" + return self.get_deriv_funcs("NoiseComponent", "toasigma") + @property_exists def d_phase_d_delay_funcs(self): """List of d_phase_d_delay functions.""" @@ -1929,7 +1943,7 @@ def d_delay_d_param_num(self, toas, param, step=1e-2): return d_delay * (u.second / unit) def d_dm_d_param(self, data, param): - """Return the derivative of dm with respect to the parameter.""" + """Return the derivative of DM with respect to the parameter.""" par = getattr(self, param) result = np.zeros(len(data)) << (u.pc / u.cm**3 / par.units) dm_df = self.dm_derivs.get(param, None) @@ -1945,6 +1959,23 @@ def d_dm_d_param(self, data, param): ) return result + def d_toasigma_d_param(self, data, param): + """Return the derivative of the scaled TOA uncertainty with respect to the parameter.""" + par = getattr(self, param) + result = np.zeros(len(data)) << (u.s / par.units) + sigma_df = self.toasigma_derivs.get(param, None) + if sigma_df is None: + if param not in self.params: # Maybe add differentiable params + raise AttributeError(f"Parameter {param} does not exist") + else: + return result + + for df in sigma_df: + result += df(data, param).to( + result.unit, equivalencies=u.dimensionless_angles() + ) + return result + def designmatrix(self, toas, acc_delay=None, incfrozen=False, incoffset=True): """Return the design matrix. diff --git a/src/pint/residuals.py b/src/pint/residuals.py index 736c5798b..ea7ad724b 100644 --- a/src/pint/residuals.py +++ b/src/pint/residuals.py @@ -500,10 +500,14 @@ def calc_time_resids( ) return (phase_resids / self.get_PSR_freq(calctype=calctype)).to(u.s) - def _calc_gls_chi2(self): + def _calc_gls_chi2(self, lognorm=False): """Compute the chi2 when correlated noise is present in the timing model. If the system is not singular, it uses Cholesky factorization to evaluate this. - If the system is singular, it uses singular value decomposition instead.""" + If the system is singular, it uses singular value decomposition instead. + + If `lognorm=True` is given, the log-normalization-factor of the likelihood + function will also be returned. + """ self.update() residuals = self.time_resids.to(u.s).value @@ -515,7 +519,11 @@ def _calc_gls_chi2(self): # For implicit offset subtraction if "PHOFF" not in self.model: M = np.append(M, np.ones((len(self.toas), 1)), axis=1) - phiinv = np.append(phiinv, [0]) + + # 1e-40 is used here instead of 0 to avoid an + # un-normalizable distribution. + # ENTERPRISE uses the same value. + phiinv = np.append(phiinv, [1e-40]) M, norm = normalize_designmatrix(M, None) @@ -529,6 +537,7 @@ def _calc_gls_chi2(self): try: c = cho_factor(mtcm) xhat = cho_solve(c, mtcy) + svd = False except LinAlgError as e: log.warning( f"Degenerate conditions encountered when computing chi-squared: {e}" @@ -541,13 +550,26 @@ def _calc_gls_chi2(self): xhat = np.dot(Vt.T, np.dot(U.T, mtcy) / s) + svd = True + newres = residuals - np.dot(M, xhat) chi2 = np.dot(newres, cinv * newres) + np.dot(xhat, phiinv * xhat) - return chi2 + if not lognorm: + return chi2 + else: + logdet_N = np.sum(np.log(np.abs(Nvec))) + logdet_Phiinv = np.sum(np.log(np.abs(phiinv))) + logdet_Sigma = ( + np.sum(np.log(np.abs(np.diag(c[0])))) + if not svd + else np.sum(np.log(np.abs(s))) + ) + log_norm = 0.5 * (logdet_N - logdet_Phiinv + logdet_Sigma) + return chi2, log_norm - def calc_chi2(self): + def calc_chi2(self, lognorm=False): """Return the weighted chi-squared for the model and toas. If the errors on the TOAs are independent this is a straightforward @@ -565,14 +587,29 @@ def calc_chi2(self): Handling of problematic results - degenerate conditions explored by a minimizer for example - may need to be checked to confirm that they correctly return infinity. + + If `lognorm=True` is given, the log-normalization-factor of the likelihood + function will also be returned. + + Parameters + ---------- + lognorm: bool + If True, return the the log-normalization-factor of the likelihood + function along with the chi2 value. + + Returns + ------- + chi2 if lognorm is False + (chi2, log_norm) if lognorm is True """ if self.model.has_correlated_errors: - return self._calc_gls_chi2() + return self._calc_gls_chi2(lognorm=lognorm) else: # Residual units are in seconds. Error units are in microseconds. toa_errors = self.get_data_error() if (toa_errors == 0.0).any(): return np.inf + # The self.time_resids is in the unit of "s", the error "us". # This is more correct way, but it is the slowest. # return (((self.time_resids / self.toas.get_errors()).decompose()**2.0).sum()).value @@ -582,10 +619,28 @@ def calc_chi2(self): # This the fastest way, but highly depend on the assumption of time_resids and # error units. Ensure only a pure number is returned. - try: - return ((self.time_resids / toa_errors.to(u.s)) ** 2.0).sum().value - except ValueError: - return ((self.time_resids / toa_errors.to(u.s)) ** 2.0).sum() + + r = self.time_resids + err = toa_errors.to(u.s) + + chi2 = ((r / err) ** 2.0).sum().value + + if not lognorm: + return chi2 + else: + log_norm = np.sum(np.log(err.value)) + return chi2, log_norm + + def lnlikelihood(self): + """Compute the log-likelihood for the model and TOAs.""" + chi2, log_norm = self.calc_chi2(lognorm=True) + return -(chi2 / 2 + log_norm) + + def d_lnlikelihood_d_whitenoise_param(self, param): + r = self.time_resids + sigma = self.get_data_error() + d_sigma_d_param = self.model.d_toasigma_d_param(self.toas, param) + return np.sum(((r / sigma) ** 2 - 1) / sigma * d_sigma_d_param) def ecorr_average(self, use_noise_model=True): """Uses the ECORR noise model time-binning to compute "epoch-averaged" residuals. diff --git a/src/pint/simulation.py b/src/pint/simulation.py index a8dcfeeca..cb64ef3b6 100644 --- a/src/pint/simulation.py +++ b/src/pint/simulation.py @@ -455,7 +455,18 @@ def make_fake_toas_fromtim( -------- :func:`make_fake_toas` """ - input_ts = pint.toa.get_TOAs(timfile, planets=model.PLANET_SHAPIRO.value) + ephem = ( + model.EPHEM.value + if hasattr(model, "EPHEM") and model.EPHEM.value is not None + else None + ) + planets = ( + model.PLANET_SHAPIRO.value + if hasattr(model, "PLANET_SHAPIRO") and model.PLANET_SHAPIRO.value is not None + else False + ) + + input_ts = pint.toa.get_TOAs(timfile, ephem=ephem, planets=planets) if input_ts.is_wideband(): dm_errors = input_ts.get_dm_errors() diff --git a/src/pint/utils.py b/src/pint/utils.py index 9fcdb31dc..a06ba5e23 100644 --- a/src/pint/utils.py +++ b/src/pint/utils.py @@ -2601,3 +2601,32 @@ def normalize_designmatrix(M, params): norm[norm == 0] = 1 return M / norm, norm + + +def akaike_information_criterion(model, toas): + """Compute the Akaike information criterion. + + Parameters + ---------- + model: pint.models.timing_model.TimingModel + The timing model + toas: pint.toas.TOAs + TOAs + + Returns + ------- + aic: float + The Akaike information criterion + """ + from pint.residuals import Residuals + + if not toas.is_wideband(): + k = ( + len(model.free_params) + if "PhaseOffset" in model.components + else len(model.free_params) + 1 + ) + lnL = Residuals(toas, model).lnlikelihood() + return 2 * (k - lnL) + else: + raise NotImplementedError diff --git a/tests/test_ddgr.py b/tests/test_ddgr.py index d436a6f7b..8925ece05 100644 --- a/tests/test_ddgr.py +++ b/tests/test_ddgr.py @@ -163,7 +163,7 @@ def test_ddgrfit_noMTOT(self): assert fDDGR.resids.calc_chi2() > chi2DDGR fDDGR.fit_toas() assert np.isclose(fDDGR.resids.calc_chi2(), chi2DDGR, atol=0.1) - assert np.isclose(fDDGR.model.M2.quantity, M2) + assert np.isclose(fDDGR.model.M2.quantity, M2, atol=1e-7) def test_ddgrfit(self): t = pint.simulation.make_fake_toas_uniform( diff --git a/tests/test_noise_models.py b/tests/test_noise_models.py index acf74afad..7165c772f 100644 --- a/tests/test_noise_models.py +++ b/tests/test_noise_models.py @@ -4,9 +4,11 @@ import numpy as np import pytest from pint.config import examplefile -from pint.models import get_model_and_toas +from pint.models import get_model_and_toas, get_model from pint.models.timing_model import Component from pint.models.noise_model import NoiseComponent +from pint.simulation import make_fake_toas_uniform +from io import StringIO noise_component_labels = [ @@ -119,3 +121,36 @@ def test_noise_basis_weights_funcs(model_and_toas, component_label): weights_ = component.get_noise_weights(toas) assert np.allclose(basis_, basis) and np.allclose(weights, weights_) + + +def test_white_noise_model_derivs(): + par = """ + ELAT 1.3 1 + ELONG 2.5 1 + F0 100 1 + F1 1e-13 1 + PEPOCH 55000 + EFAC mjd 49999 53000 2 1 + EQUAD mjd 49999 53000 0.5 1 + EFAC mjd 53000 55000 1.5 1 + EQUAD mjd 53000 55000 1 1 + """ + m = get_model(StringIO(par)) + t = make_fake_toas_uniform(50000, 55000, 1000, m, add_noise=True) + + dq1 = m.d_toasigma_d_param(t, "EQUAD1") + dq2 = m.d_toasigma_d_param(t, "EQUAD2") + df1 = m.d_toasigma_d_param(t, "EFAC1") + df2 = m.d_toasigma_d_param(t, "EFAC2") + + sigma = m.scaled_toa_uncertainty(t) + + assert np.isclose(df1[0], sigma[0] / m.EFAC1.quantity) + assert np.isclose(df2[-1], sigma[-1] / m.EFAC2.quantity) + + assert np.isclose( + dq1[0], sigma[0] * m.EQUAD1.quantity * (m.EFAC1.quantity / sigma[0]) ** 2 + ) + assert np.isclose( + dq2[-1], sigma[-1] * m.EQUAD2.quantity * (m.EFAC2.quantity / sigma[-1]) ** 2 + ) diff --git a/tests/test_noisefit.py b/tests/test_noisefit.py new file mode 100644 index 000000000..b2690d7c2 --- /dev/null +++ b/tests/test_noisefit.py @@ -0,0 +1,61 @@ +from pint.models import get_model +from pint.simulation import make_fake_toas_uniform +from pint.fitter import DownhillWLSFitter + +from io import StringIO +import numpy as np + +par = """ + ELAT 1.3 1 + ELONG 2.5 1 + F0 100 1 + F1 1e-13 1 + PEPOCH 55000 + EPHEM DE440 + EFAC mjd 50000 53000 2 1 + EQUAD mjd 53000 55000 0.8 1 +""" + +m = get_model(StringIO(par)) +t = make_fake_toas_uniform(50000, 55000, 200, m, add_noise=True) + + +def test_white_noise_fit(): + assert m.EFAC1.uncertainty_value == 0 and m.EQUAD1.uncertainty_value == 0 + + ftr = DownhillWLSFitter(t, m) + ftr.fit_toas(maxiter=5) + + assert ( + ftr.model.EFAC1.uncertainty_value > 0 and ftr.model.EQUAD1.uncertainty_value > 0 + ) + assert ( + np.abs(m.EFAC1.value - ftr.model.EFAC1.value) + / ftr.model.EFAC1.uncertainty_value + < 3 + ) + assert ( + np.abs(m.EQUAD1.value - ftr.model.EQUAD1.value) + / ftr.model.EQUAD1.uncertainty_value + < 3 + ) + + +def test_white_noise_refit(): + ftr = DownhillWLSFitter(t, m) + + ftr.model.EFAC1.value = 1.5 + ftr.model.EQUAD1.value = 0.5 + + ftr.fit_toas(maxiter=5) + + assert ( + np.abs(m.EFAC1.value - ftr.model.EFAC1.value) + / ftr.model.EFAC1.uncertainty_value + < 3 + ) + assert ( + np.abs(m.EQUAD1.value - ftr.model.EQUAD1.value) + / ftr.model.EQUAD1.uncertainty_value + < 3 + ) diff --git a/tests/test_utils.py b/tests/test_utils.py index 18988b4d0..b2e653757 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -56,6 +56,7 @@ print_color_examples, parse_time, info_string, + akaike_information_criterion, ) @@ -871,3 +872,10 @@ def test_parse_time(t): def test_info_str(): info = info_string() dinfo = info_string(detailed=True) + + +def test_aic(): + m = tm.get_model(os.path.join(datadir, "B1855+09_NANOGrav_9yv1.gls.par")) + t = toa.get_TOAs(os.path.join(datadir, "B1855+09_NANOGrav_9yv1.tim")) + + assert np.isfinite(akaike_information_criterion(m, t)) diff --git a/tox.ini b/tox.ini index f6d4241a0..05ed53ca3 100644 --- a/tox.ini +++ b/tox.ini @@ -91,6 +91,7 @@ deps = sphinx_rtd_theme pygments jupyter + numdifftools nbconvert pytest jupytext