diff --git a/hierarc/Likelihood/LensLikelihood/td_mag_likelihood.py b/hierarc/Likelihood/LensLikelihood/td_mag_likelihood.py index 750b0de..13c0388 100644 --- a/hierarc/Likelihood/LensLikelihood/td_mag_likelihood.py +++ b/hierarc/Likelihood/LensLikelihood/td_mag_likelihood.py @@ -7,6 +7,8 @@ class TDMagLikelihood(object): """ likelihood of time delays and magnification likelihood + This likelihood uses linear flux units and linear lensing magnifications. + """ def __init__(self, time_delay_measured, cov_td_measured, amp_measured, cov_amp_measured, fermat_diff, magnification_model, cov_model, magnitude_zero_point=20): diff --git a/hierarc/Likelihood/LensLikelihood/td_mag_magnitude_likelihood.py b/hierarc/Likelihood/LensLikelihood/td_mag_magnitude_likelihood.py new file mode 100644 index 0000000..8555193 --- /dev/null +++ b/hierarc/Likelihood/LensLikelihood/td_mag_magnitude_likelihood.py @@ -0,0 +1,94 @@ +import numpy as np +from lenstronomy.Util import constants as const +from lenstronomy.Util.data_util import magnitude2cps + + +class TDMagMagnitudeLikelihood(object): + """ + likelihood of time delays and magnification likelihood + + This likelihood uses astronomical magnitude units in flux measurement and lensing magnification + and Gaussian uncertainties in this space. + + """ + + def __init__(self, time_delay_measured, cov_td_measured, magnitude_measured, cov_magnitude_measured, + fermat_diff, magnification_model, cov_model): + """ + + :param time_delay_measured: array, relative time delays (relative to the first image) [days] + :param cov_td_measured: 2d array, error covariance matrix of time delay measurement [days^2] + :param magnitude_measured: array, astronomical magnitude of measured fluxes of image positions + :param cov_magnitude_measured: 2d array, error covariance matrix of the measured amplitudes + :param fermat_diff: mean Fermat potential differences (relative to the first image) in arcsec^2 + :param magnification_model: mean lensing magnification of the model prediction + in units of astronomical magnitudes (array of length of the images, (mean of - 2.5 * log10(mu)) + :param cov_model: 2d array (length relative time delays + image magnifications); + model fermat potential differences and lensing magnification in astronomical magnitudes covariances + """ + + self._data_vector = np.append(time_delay_measured, magnitude_measured) + self._cov_td_measured = np.array(cov_td_measured) + self._cov_magnitude_measured = np.array(cov_magnitude_measured) + # check sizes of covariances matches + n_tot = len(self._data_vector) + self._n_td = len(time_delay_measured) + self._n_amp = len(magnitude_measured) + assert self._n_td == len(cov_td_measured) + assert self._n_amp == len(cov_magnitude_measured) + assert n_tot == len(cov_model) + # merge data covariance matrices from time delay and image amplitudes + self._cov_data = np.zeros((n_tot, n_tot)) + self._cov_data[:self._n_td, :self._n_td] = self._cov_td_measured + self._cov_data[self._n_td:, self._n_td:] = self._cov_magnitude_measured + # self._fermat_diff = fermat_diff # in units arcsec^2 + self._fermat_unit_conversion = const.Mpc / const.c / const.day_s * const.arcsec ** 2 + # self._mag_model = mag_model + self._model_tot = np.append(fermat_diff, magnification_model) + self._cov_model = cov_model + self.num_data = n_tot + + def log_likelihood(self, ddt, mu_intrinsic): + """ + + :param ddt: time-delay distance (physical Mpc) + :param mu_intrinsic: intrinsic brightness of the source (already incorporating the inverse MST transform) + :return: log likelihood of the measured magnified images given the source brightness + """ + model_vector, cov_tot = self._model_cov(ddt, mu_intrinsic) + # invert matrix + try: + cov_tot_inv = np.linalg.inv(cov_tot) + except: + return -np.inf + # difference to data vector + delta = self._data_vector - model_vector + # evaluate likelihood + lnlikelihood = -delta.dot(cov_tot_inv.dot(delta)) / 2. + sign_det, lndet = np.linalg.slogdet(cov_tot) + lnlikelihood -= 1 / 2. * (self.num_data * np.log(2 * np.pi) + lndet) + return lnlikelihood + + def _model_cov(self, ddt, mu_intrinsic): + """ + combined covariance matrix of the data and model when marginialized over the Gaussian model uncertainties + in the Fermat potential and magnification. + + :param ddt: time-delay distance (physical Mpc) + :param mu_intrinsic: intrinsic brightness of the source (already incorporating the inverse MST transform) + :return: model vector, combined covariance matrix + """ + # compute model predicted magnified image amplitude and time delay + + model_scale = np.append(ddt * self._fermat_unit_conversion * np.ones(self._n_td), np.ones(self._n_amp)) + model_vector = np.zeros_like(self._model_tot) + # time delay prediction + model_vector[:self._n_td] = ddt * self._fermat_unit_conversion * self._model_tot[:self._n_td] + # lensed astronomical magnitude prediction + model_vector[self._n_td:] = self._model_tot[self._n_td:] + mu_intrinsic + # scale model covariance matrix with model_scale vector (in quadrature), + # shift in astronomical magnitudes do not change covariance matrix + cov_model = model_scale * (self._cov_model * model_scale).T + # combine data and model covariance matrix + cov_tot = self._cov_data + cov_model + return model_vector, cov_tot diff --git a/test/test_Likelihood/test_LensLikelihood/test_td_mag_magnitude_likelihood.py b/test/test_Likelihood/test_LensLikelihood/test_td_mag_magnitude_likelihood.py new file mode 100644 index 0000000..9003cc4 --- /dev/null +++ b/test/test_Likelihood/test_LensLikelihood/test_td_mag_magnitude_likelihood.py @@ -0,0 +1,74 @@ +import pytest +import numpy as np +import numpy.testing as npt +from lenstronomy.Util import constants as const +from hierarc.Likelihood.LensLikelihood.td_mag_likelihood import TDMagLikelihood +from hierarc.Likelihood.LensLikelihood.td_mag_magnitude_likelihood import TDMagMagnitudeLikelihood +from lenstronomy.Util.data_util import magnitude2cps + + +class TestMagnificationLikelihood(object): + + def setup(self): + pass + + def test_log_likelihood(self): + ddt = 1000 # time-delay distance in units Mpc + num = 4 + magnification_model = np.ones(num) * 4 + magnitude_intrinsic = 19 + + magnitude_zero_point = 20 + amp_int = magnitude2cps(magnitude=magnitude_intrinsic, magnitude_zero_point=magnitude_zero_point) + + amp_measured = magnification_model * amp_int + rel_sigma = 0.1 + cov_amp_measured = np.diag((amp_measured*rel_sigma)**2) + magnification_model_magnitude = - 2.5 * np.log10(magnification_model) + magnitude_measured = magnitude_intrinsic - 2.5 * np.log10(magnification_model) + cov_magnitude_measured = np.diag(rel_sigma * 1.086 * np.ones(num)) # translating relative scatter in linear flux units to astronomical magnitudes + + time_delay_measured = np.ones(num - 1) * 10 + cov_td_measured = np.ones((num-1, num-1)) + fermat_unit_conversion = const.Mpc / const.c / const.day_s * const.arcsec ** 2 + fermat_diff = time_delay_measured / fermat_unit_conversion / ddt + model_vector = np.append(fermat_diff, magnification_model) + cov_model = np.diag((model_vector/10)**2) + + cov_model_magnitude = np.diag(np.append(fermat_diff * rel_sigma, rel_sigma * 1.086 * np.ones(num))) + # un-normalized likelihood, comparing linear flux likelihood and magnification likelihood + + # linear flux likelihood + likelihood = TDMagLikelihood(time_delay_measured, cov_td_measured, amp_measured, cov_amp_measured, + fermat_diff, magnification_model, cov_model, magnitude_zero_point=magnitude_zero_point) + + logl = likelihood.log_likelihood(ddt=ddt, mu_intrinsic=magnitude_intrinsic) + model_vector, cov_tot = likelihood._model_cov(ddt, mu_intrinsic=magnitude_intrinsic) + sign_det, lndet = np.linalg.slogdet(cov_tot) + logl_norm = -1 / 2. * (likelihood.num_data * np.log(2 * np.pi) + lndet) + npt.assert_almost_equal(logl, logl_norm, decimal=6) + + # astronomical magnitude likelihood + likelihood_magnitude = TDMagMagnitudeLikelihood(time_delay_measured, cov_td_measured, magnitude_measured, + cov_magnitude_measured, fermat_diff, magnification_model_magnitude, + cov_model_magnitude) + logl_magnitude = likelihood_magnitude.log_likelihood(ddt=ddt, mu_intrinsic=magnitude_intrinsic) + npt.assert_almost_equal(logl_magnitude - logl, 0, decimal=1) + + num = 4 + magnification_model = np.ones(num) + cov_td_measured = np.zeros((num - 1, num - 1)) + cov_magnitude_measured = np.zeros((num, num)) + amp_int = 10 + magnitude_measured = magnitude_intrinsic - 2.5 * np.log10(magnification_model) + cov_model = np.zeros((num + num - 1, num + num - 1)) + + likelihood = TDMagMagnitudeLikelihood(time_delay_measured, cov_td_measured, magnitude_measured, + cov_magnitude_measured, fermat_diff, magnification_model_magnitude, + cov_model) + logl = likelihood.log_likelihood(ddt=ddt, mu_intrinsic=1) + assert logl == -np.inf + + +if __name__ == '__main__': + pytest.main()