diff --git a/AUTHORS.rst b/AUTHORS.rst index 278659f69..3112065c8 100644 --- a/AUTHORS.rst +++ b/AUTHORS.rst @@ -9,7 +9,7 @@ Current maintainers * Anowar Shajib `ajshajib `_ * Daniel Gilman `dangilman `_ -lenstronomy developer `mailing list `_ +Contact the lenstronomy developers via `email `_ if you have questions. @@ -28,6 +28,7 @@ Contributors (alphabetic) * Aymeric Galan `aymgal `_ * Matthew R. Gomer `mattgomer `_ * Natalie B. Hogg `nataliehogg `_ +* Tyler Hughes * Daniel Johnson `DanJohnson98 `_ * Felix A. Kuhn * Felix Mayor @@ -37,6 +38,7 @@ Contributors (alphabetic) * Brian Nord `bnord `_ * Giulia Pagano * Ji Won Park `jiwoncpark `_ +* Jackson O'Donnell `jhod0 `_ * Thomas Schmidt `Thomas-01 `_ * Dominique Sluse * Luca Teodori `lucateo `_ @@ -53,5 +55,12 @@ Contributors (alphabetic) Past development lead --------------------- -The initial source code of lenstronomy was developed by Simon Birrer `sibirrer `_ -in 2014-2018 and made public in 2018. The lenstronomy developement moved to the project repository in 2022. \ No newline at end of file +The initial source code of lenstronomy was developed by Simon Birrer (`sibirrer `_) +in 2014-2018 and made public in 2018. From 2018-2022 the development of lenstronomy was hosted on Simon Birrer's +repository with increased contributions from many people. +The lenstronomy development moved to the `project repository `_ in 2022. + + +Lenstronomy logo +---------------- +The lenstronomy logo was designed by Zoe Alexander `zoe-blyss `_. diff --git a/PUBLISHED.rst b/PUBLISHED.rst index 60882370d..5efaa72d3 100644 --- a/PUBLISHED.rst +++ b/PUBLISHED.rst @@ -2,10 +2,10 @@ Published work with lenstronomy =============================== -In this section you can find the concept papers **lenstronomy** is based on the list of science publications that made -use of **lenstronomy**. Please let the developers know when you publish a paper that made use of **lenstronomy**. -We are happy to include your publication in this list. - +In this section you can find the concept papers **lenstronomy** is based on and a list of science publications that made +use of **lenstronomy** before 09/2022. +For a more complete and current list of publications using lenstronomy we refer to the `NASA/ADS query `_ +(this incudes all publications citing lenstronomy papers, which is not the same as publications making active use of the software). Core lenstronomy methodology and software publications @@ -44,6 +44,10 @@ Related software publications +===================================== +Scientific publication before 09/2022 +===================================== + @@ -235,6 +239,10 @@ Galaxy formation and evolution * Early results from GLASS-JWST. V: the first rest-frame optical size-luminosity relation of galaxies at z>7; `Yang et al. 2022b `_ *galaxy size measurement from JWST data with Galight/lenstronomy* +* A New Polar Ring Galaxy Discovered in the COSMOS Field; `Nishimura et al. 2022 `_ + +* Webb's PEARLS: dust attenuation and gravitational lensing in the backlit-galaxy system VV 191; `Keel et al. 2022 `_ + Automatized Lens Modeling @@ -347,7 +355,6 @@ Large scale structure - Others ------ diff --git a/README.rst b/README.rst index 19618f3a2..35266a40c 100644 --- a/README.rst +++ b/README.rst @@ -1,7 +1,6 @@ -==================================================== -lenstronomy - gravitational lensing software package -==================================================== +.. image:: https://raw.githubusercontent.com/lenstronomy/lenstronomy/main/docs/figures/logo_text.png + :target: https://raw.githubusercontent.com/lenstronomy/lenstronomy/main/docs/figures/logo_text.png .. image:: https://github.com/lenstronomy/lenstronomy/workflows/Tests/badge.svg :target: https://github.com/lenstronomy/lenstronomy/actions @@ -26,8 +25,9 @@ lenstronomy - gravitational lensing software package .. image:: https://img.shields.io/badge/arXiv-1803.09746%20-yellowgreen.svg :target: https://arxiv.org/abs/1803.09746 -.. image:: https://raw.githubusercontent.com/lenstronomy/lenstronomy/main/docs/figures/readme_fig.png - :target: https://raw.githubusercontent.com/lenstronomy/lenstronomy/main/docs/figures/readme_fig.png +.. + .. image:: https://raw.githubusercontent.com/lenstronomy/lenstronomy/main/docs/figures/readme_fig.png + :target: https://raw.githubusercontent.com/lenstronomy/lenstronomy/main/docs/figures/readme_fig.png ``lenstronomy`` is a multi-purpose software package to model strong gravitational lenses. diff --git a/docs/conf.py b/docs/conf.py index 402c636d0..4b38aae57 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -31,7 +31,13 @@ # Add any Sphinx extension module names here, as strings. They can be extensions # coming with Sphinx (named 'sphinx.ext.*') or your custom ones. -extensions = ['sphinx.ext.autodoc', 'sphinx.ext.coverage', 'sphinx.ext.mathjax', 'sphinx.ext.viewcode'] +extensions = ['sphinx.ext.autodoc', 'sphinx.ext.coverage', 'sphinx.ext.mathjax', 'sphinx.ext.viewcode']\ +# , 'sphinx.ext.autosectionlabel'] + +autodoc_default_options = { + 'member-order': 'bysource', + 'special-members': '__init__', +} # Add any paths that contain templates here, relative to this directory. templates_path = ['_templates'] diff --git a/docs/figures/logo_text.png b/docs/figures/logo_text.png new file mode 100644 index 000000000..f4a828d9c Binary files /dev/null and b/docs/figures/logo_text.png differ diff --git a/docs/index.rst b/docs/index.rst index 593c8e48c..0fd0fada6 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -12,11 +12,11 @@ Contents: :maxdepth: 2 installation - modules usage contributing mailinglist authors published affiliatedpackages + lenstronomy history diff --git a/docs/lenstronomy.Data.rst b/docs/lenstronomy.Data.rst index aadffa78e..b2e47e67c 100644 --- a/docs/lenstronomy.Data.rst +++ b/docs/lenstronomy.Data.rst @@ -1,38 +1,53 @@ -lenstronomy\.Data package -========================= +lenstronomy.Data package +======================== Submodules ---------- -lenstronomy\.Data\.coord\_transforms module -------------------------------------------- +lenstronomy.Data.coord\_transforms module +----------------------------------------- .. automodule:: lenstronomy.Data.coord_transforms - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: -lenstronomy\.Data\.imaging\_data module ---------------------------------------- +lenstronomy.Data.image\_noise module +------------------------------------ + +.. automodule:: lenstronomy.Data.image_noise + :members: + :undoc-members: + :show-inheritance: + +lenstronomy.Data.imaging\_data module +------------------------------------- .. automodule:: lenstronomy.Data.imaging_data - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: -lenstronomy\.Data\.psf module ------------------------------ +lenstronomy.Data.pixel\_grid module +----------------------------------- -.. automodule:: lenstronomy.Data.psf - :members: - :undoc-members: - :show-inheritance: +.. automodule:: lenstronomy.Data.pixel_grid + :members: + :undoc-members: + :show-inheritance: +lenstronomy.Data.psf module +--------------------------- + +.. automodule:: lenstronomy.Data.psf + :members: + :undoc-members: + :show-inheritance: Module contents --------------- .. automodule:: lenstronomy.Data - :members: - :undoc-members: - :show-inheritance: + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/lenstronomy.Sampling.rst b/docs/lenstronomy.Sampling.rst index 1aeaacc8c..08a211d8a 100644 --- a/docs/lenstronomy.Sampling.rst +++ b/docs/lenstronomy.Sampling.rst @@ -45,6 +45,14 @@ lenstronomy.Sampling.special\_param module :undoc-members: :show-inheritance: +lenstronomy.Sampling.param\_group module +------------------------------------------ + +.. automodule:: lenstronomy.Sampling.param_group + :members: + :undoc-members: + :show-inheritance: + Module contents --------------- diff --git a/docs/usage.rst b/docs/usage.rst index c33f7b182..509610d4e 100644 --- a/docs/usage.rst +++ b/docs/usage.rst @@ -17,7 +17,7 @@ user to directly access a lot of tools and each module can also be used as stand Example notebooks ----------------- -We have made an extension module available at `https://github.com/lenstronomy/lenstronom-tutorials `_. +We have made an extension module available at `https://github.com/lenstronomy/lenstronomy-tutorials `_. You can find simple example notebooks for various cases. The latest versions of the notebooks should be compatible with the recent pip version of lenstronomy. diff --git a/lenstronomy/Analysis/kinematics_api.py b/lenstronomy/Analysis/kinematics_api.py index 52b1fdd71..8e03c1ef0 100644 --- a/lenstronomy/Analysis/kinematics_api.py +++ b/lenstronomy/Analysis/kinematics_api.py @@ -51,7 +51,7 @@ def __init__(self, z_lens, z_source, kwargs_model, kwargs_aperture, kwargs_seein :param kwargs_mge_mass: keyword arguments that go into the MGE decomposition routine :param kwargs_mge_light: keyword arguments that go into the MGE decomposition routine :param sampling_number: int, number of spectral rendering to compute the light weighted integrated LOS - dispersion within the aperture. This keyword should be chosen high enough to result in converged results within the tolerance. + dispersion within the aperture. This keyword should be chosen high enough to result in converged results within the tolerance. :param num_kin_sampling: number of kinematic renderings on a total IFU :param num_psf_sampling: number of PSF displacements for each kinematic rendering on the IFU """ diff --git a/lenstronomy/Cosmo/background.py b/lenstronomy/Cosmo/background.py index b05411483..30e5df32a 100644 --- a/lenstronomy/Cosmo/background.py +++ b/lenstronomy/Cosmo/background.py @@ -17,7 +17,7 @@ def __init__(self, cosmo=None, interp=False, **kwargs_interp): :param cosmo: instance of astropy.cosmology :param interp: boolean, if True, uses interpolated cosmology to evaluate specific redshifts :param kwargs_interp: keyword arguments of CosmoInterp specifying the interpolation interval and maximum - redshift + redshift :return: Background class with instance of astropy.cosmology """ @@ -29,9 +29,11 @@ def __init__(self, cosmo=None, interp=False, **kwargs_interp): else: self.cosmo = cosmo - def a_z(self, z): + @staticmethod + def a_z(z): """ returns scale factor (a_0 = 1) for given redshift + :param z: redshift :return: scale factor """ @@ -72,6 +74,7 @@ def T_xy(self, z_observer, z_source): def rho_crit(self): """ critical density + :return: value in M_sol/Mpc^3 """ h = self.cosmo.H(0).value / 100. diff --git a/lenstronomy/Cosmo/cosmo_solver.py b/lenstronomy/Cosmo/cosmo_solver.py index 339e7db18..e7e6b8cdb 100644 --- a/lenstronomy/Cosmo/cosmo_solver.py +++ b/lenstronomy/Cosmo/cosmo_solver.py @@ -85,15 +85,20 @@ class InvertCosmo(object): """ class to do an interpolation and call the inverse of this interpolation to get H_0 and omega_m """ - def __init__(self, z_d, z_s, H0_range=np.linspace(10, 100, 90), omega_m_range=np.linspace(0.05, 1, 95)): + def __init__(self, z_d, z_s, H0_range=None, omega_m_range=None): self.z_d = z_d self.z_s = z_s + if H0_range is None: + H0_range = np.linspace(10, 100, 90) + if omega_m_range is None: + omega_m_range = np.linspace(0.05, 1, 95) self._H0_range = H0_range self._omega_m_range = omega_m_range def _make_interpolation(self): """ creates an interpolation grid in H_0, omega_m and computes quantities in Dd and Ds_Dds + :return: """ grid2d = np.dstack(np.meshgrid(self._H0_range, self._omega_m_range)).reshape(-1, 2) @@ -115,6 +120,7 @@ def _make_interpolation(self): def get_cosmo(self, Dd, Ds_Dds): """ return the values of H0 and omega_m computed with an interpolation + :param Dd: flat :param Ds_Dds: float :return: diff --git a/lenstronomy/Cosmo/kde_likelihood.py b/lenstronomy/Cosmo/kde_likelihood.py index 1d4dea836..dbed4d4d0 100644 --- a/lenstronomy/Cosmo/kde_likelihood.py +++ b/lenstronomy/Cosmo/kde_likelihood.py @@ -14,11 +14,11 @@ def __init__(self, D_d_sample, D_delta_t_sample, kde_type='scipy_gaussian', band :param D_d_sample: 1-d numpy array of angular diameter distances to the lens plane :param D_delta_t_sample: 1-d numpy array of time-delay distances - kde_type : string - The kernel to use. Valid kernels are - 'scipy_gaussian' or - ['gaussian'|'tophat'|'epanechnikov'|'exponential'|'linear'|'cosine'] - Default is 'gaussian'. + :param kde_type: The kernel to use. Valid kernels are + 'scipy_gaussian' or + ['gaussian'|'tophat'|'epanechnikov'|'exponential'|'linear'|'cosine'] + Default is 'gaussian'. + :type kde_type: string :param bandwidth: width of kernel (in same units as the angular diameter quantities) """ values = np.vstack([D_d_sample, D_delta_t_sample]) diff --git a/lenstronomy/Cosmo/lcdm.py b/lenstronomy/Cosmo/lcdm.py index 491313652..618f5d21e 100644 --- a/lenstronomy/Cosmo/lcdm.py +++ b/lenstronomy/Cosmo/lcdm.py @@ -40,6 +40,7 @@ def _get_cosom(self, H_0, Om0, Ode0=None): def D_d(self, H_0, Om0, Ode0=None): """ angular diameter to deflector + :param H_0: Hubble parameter [km/s/Mpc] :param Om0: normalized matter density at present time :return: float [Mpc] @@ -50,6 +51,7 @@ def D_d(self, H_0, Om0, Ode0=None): def D_s(self, H_0, Om0, Ode0=None): """ angular diameter to source + :param H_0: Hubble parameter [km/s/Mpc] :param Om0: normalized matter density at present time :return: float [Mpc] @@ -60,6 +62,7 @@ def D_s(self, H_0, Om0, Ode0=None): def D_ds(self, H_0, Om0, Ode0=None): """ angular diameter from deflector to source + :param H_0: Hubble parameter [km/s/Mpc] :param Om0: normalized matter density at present time :return: float [Mpc] @@ -69,7 +72,8 @@ def D_ds(self, H_0, Om0, Ode0=None): def D_dt(self, H_0, Om0, Ode0=None): """ - time delay distance + time-delay distance + :param H_0: Hubble parameter [km/s/Mpc] :param Om0: normalized matter density at present time :return: float [Mpc] diff --git a/lenstronomy/Cosmo/lens_cosmo.py b/lenstronomy/Cosmo/lens_cosmo.py index ee314ef9d..93af64b4c 100644 --- a/lenstronomy/Cosmo/lens_cosmo.py +++ b/lenstronomy/Cosmo/lens_cosmo.py @@ -30,6 +30,7 @@ def __init__(self, z_lens, z_source, cosmo=None): def a_z(self, z): """ convert redshift into scale factor + :param z: redshift :return: scale factor """ @@ -75,6 +76,7 @@ def ddt(self): def sigma_crit(self): """ returns the critical projected lensing mass density in units of M_sun/Mpc^2 + :return: critical projected lensing mass density """ if not hasattr(self, '_sigma_crit_mpc'): @@ -89,6 +91,7 @@ def sigma_crit_angle(self): """ returns the critical surface density in units of M_sun/arcsec^2 (in physical solar mass units) when provided a physical mass per physical Mpc^2 + :return: critical projected mass density """ if not hasattr(self, '_sigma_crit_arcsec'): @@ -101,6 +104,7 @@ def sigma_crit_angle(self): def phys2arcsec_lens(self, phys): """ convert physical Mpc into arc seconds + :param phys: physical distance [Mpc] :return: angular diameter [arcsec] """ @@ -109,6 +113,7 @@ def phys2arcsec_lens(self, phys): def arcsec2phys_lens(self, arcsec): """ convert angular to physical quantities for lens plane + :param arcsec: angular size at lens plane [arcsec] :return: physical size at lens plane [Mpc] """ @@ -117,6 +122,7 @@ def arcsec2phys_lens(self, arcsec): def arcsec2phys_source(self, arcsec): """ convert angular to physical quantities for source plane + :param arcsec: angular size at source plane [arcsec] :return: physical size at source plane [Mpc] """ @@ -125,6 +131,7 @@ def arcsec2phys_source(self, arcsec): def kappa2proj_mass(self, kappa): """ convert convergence to projected mass M_sun/Mpc^2 + :param kappa: lensing convergence :return: projected mass [M_sun/Mpc^2] """ @@ -133,6 +140,7 @@ def kappa2proj_mass(self, kappa): def mass_in_theta_E(self, theta_E): """ mass within Einstein radius (area * epsilon crit) [M_sun] + :param theta_E: Einstein radius [arcsec] :return: mass within Einstein radius [M_sun] """ @@ -225,6 +233,7 @@ def nfw_M_theta_r200(self, M): def sis_theta_E2sigma_v(self, theta_E): """ converts the lensing Einstein radius into a physical velocity dispersion + :param theta_E: Einstein radius (in arcsec) :return: velocity dispersion in units (km/s) """ @@ -234,6 +243,7 @@ def sis_theta_E2sigma_v(self, theta_E): def sis_sigma_v2theta_E(self, v_sigma): """ converts the velocity dispersion into an Einstein radius for a SIS profile + :param v_sigma: velocity dispersion (km/s) :return: theta_E (arcsec) """ @@ -245,6 +255,7 @@ def uldm_angular2phys(self, kappa_0, theta_c): converts the anguar parameters entering the LensModel Uldm() (Ultra Light Dark Matter) class in physical masses, i.e. the total soliton mass and the mass of the particle + :param kappa_0: central convergence of profile :param theta_c: core radius (in arcseconds) :return: m_eV_log10, M_sol_log10, the log10 of the masses, m in eV and M in M_sun @@ -261,6 +272,7 @@ def uldm_mphys2angular(self, m_log10, M_log10): """ converts physical ULDM mass in the ones, in angular units, that enter the LensModel Uldm() class + :param m_log10: exponent of ULDM mass in eV :param M_log10: exponent of soliton mass in M_sun :return: kappa_0, theta_c, the central convergence and core radius (in arcseconds) @@ -275,3 +287,37 @@ def uldm_mphys2angular(self, m_log10, M_log10): theta_c = r_c / D_Lens / const.arcsec return kappa_0, theta_c + def sersic_m_star2k_eff(self, m_star, R_sersic, n_sersic): + """ + translates a total stellar mass into 'k_eff', the convergence at + 'R_sersic' (effective radius or half-light radius) for a Sersic profile + + :param m_star: total stellar mass in physical Msun + :param R_sersic: half-light radius in arc seconds + :param n_sersic: Sersic index + :return: k_eff + """ + # compute mass integral + from lenstronomy.LensModel.Profiles.sersic_utils import SersicUtil + sersic_util = SersicUtil() + norm_integral = sersic_util.total_flux(amp=1, R_sersic=R_sersic, n_sersic=n_sersic) + # compute total kappa normalization and re + k_eff = m_star / self.sigma_crit_angle + # renormalize + k_eff /= norm_integral + return k_eff + + def sersic_k_eff2m_star(self, k_eff, R_sersic, n_sersic): + """ + translates convergence at half-light radius to total integrated physical stellar mass for a Sersic profile + + :param k_eff: lensing convergence at half-light radius + :param R_sersic: half-light radius in arc seconds + :param n_sersic: Sersic index + :return: stellar mass in physical Msun + """ + from lenstronomy.LensModel.Profiles.sersic_utils import SersicUtil + sersic_util = SersicUtil() + norm_integral = sersic_util.total_flux(amp=1, R_sersic=R_sersic, n_sersic=n_sersic) + m_star = k_eff *self.sigma_crit_angle * norm_integral + return m_star diff --git a/lenstronomy/Cosmo/nfw_param.py b/lenstronomy/Cosmo/nfw_param.py index 592c6a952..687b11e75 100644 --- a/lenstronomy/Cosmo/nfw_param.py +++ b/lenstronomy/Cosmo/nfw_param.py @@ -81,6 +81,7 @@ def rho0_c(self, c, z): def c_rho0(self, rho0, z): """ computes the concentration given density normalization rho_0 in h^2/Mpc^3 (physical) (inverse of function rho0_c) + :param rho0: density normalization in h^2/Mpc^3 (physical) :param z: redshift :return: concentration parameter c diff --git a/lenstronomy/Data/coord_transforms.py b/lenstronomy/Data/coord_transforms.py index fc73e5a5d..a6f2a14d1 100644 --- a/lenstronomy/Data/coord_transforms.py +++ b/lenstronomy/Data/coord_transforms.py @@ -14,6 +14,7 @@ class to handle linear coordinate transformations of a square pixel image def __init__(self, transform_pix2angle, ra_at_xy_0, dec_at_xy_0): """ initialize the coordinate-to-pixel transform and their inverse + :param transform_pix2angle: 2x2 matrix, mapping of pixel to coordinate :param ra_at_xy_0: ra coordinate at pixel (0,0) :param dec_at_xy_0: dec coordinate at pixel (0,0) @@ -82,6 +83,7 @@ def map_pix2coord(self, x, y): def pixel_area(self): """ angular area of a pixel in the image + :return: area [arcsec^2] """ return np.abs(linalg.det(self._Mpix2a)) @@ -90,6 +92,7 @@ def pixel_area(self): def pixel_width(self): """ size of pixel + :return: sqrt(pixel_area) """ return np.sqrt(self.pixel_area) @@ -110,6 +113,7 @@ def coordinate_grid(self, nx, ny): def shift_coordinate_system(self, x_shift, y_shift, pixel_unit=False): """ shifts the coordinate system + :param x_shift: shift in x (or RA) :param y_shift: shift in y (or DEC) :param pixel_unit: bool, if True, units of pixels in input, otherwise RA/DEC @@ -121,6 +125,7 @@ def _shift_coordinates(self, x_shift, y_shift, pixel_unit=False): """ shifts the coordinate system + :param x_shift: shift in x (or RA) :param y_shift: shift in y (or DEC) :param pixel_unit: bool, if True, units of pixels in input, otherwise RA/DEC diff --git a/lenstronomy/Data/image_noise.py b/lenstronomy/Data/image_noise.py index d49535d36..138099fe7 100644 --- a/lenstronomy/Data/image_noise.py +++ b/lenstronomy/Data/image_noise.py @@ -18,10 +18,10 @@ def __init__(self, image_data, exposure_time=None, background_rms=None, noise_ma :param image_data: numpy array, pixel data values :param exposure_time: int or array of size the data; exposure time - (common for all pixels or individually for each individual pixel) + (common for all pixels or individually for each individual pixel) :param background_rms: root-mean-square value of Gaussian background noise :param noise_map: int or array of size the data; joint noise sqrt(variance) of each individual pixel. - Overwrites meaning of background_rms and exposure_time. + Overwrites meaning of background_rms and exposure_time. :param gradient_boost_factor: None or float, variance terms added in quadrature scaling with gradient^2 * gradient_boost_factor """ diff --git a/lenstronomy/Data/imaging_data.py b/lenstronomy/Data/imaging_data.py index 778047782..d86ddb750 100644 --- a/lenstronomy/Data/imaging_data.py +++ b/lenstronomy/Data/imaging_data.py @@ -29,8 +29,7 @@ class to handle the data, coordinate system and masking, including convolution w If this keyword is set, the other noise properties will be ignored. - Notes: - ------ + ** notes ** the likelihood for the data given model P(data|model) is defined in the function below. Please make sure that your definitions and units of 'exposure_map', 'background_rms' and 'image_data' are in accordance with the likelihood function. In particular, make sure that the Poisson noise contribution is defined in the count rate. @@ -43,7 +42,7 @@ def __init__(self, image_data, exposure_time=None, background_rms=None, noise_ma :param image_data: 2d numpy array of the image data :param exposure_time: int or array of size the data; exposure time - (common for all pixels or individually for each individual pixel) + (common for all pixels or individually for each individual pixel) :param background_rms: root-mean-square value of Gaussian background noise in units counts per second :param noise_map: int or array of size the data; joint noise sqrt(variance) of each individual pixel. :param gradient_boost_factor: None or float, variance terms added in quadrature scaling with diff --git a/lenstronomy/Data/psf.py b/lenstronomy/Data/psf.py index 8b7be2eb7..d6ff24a54 100644 --- a/lenstronomy/Data/psf.py +++ b/lenstronomy/Data/psf.py @@ -24,9 +24,9 @@ def __init__(self, psf_type='NONE', fwhm=None, truncation=5, pixel_size=None, ke :param kernel_point_source: 2d numpy array, odd length, centered PSF of a point source (if not normalized, will be normalized) :param psf_error_map: uncertainty in the PSF model per pixel (size of data, not super-sampled). 2d numpy array. - Size can be larger or smaller than the pixel-sized PSF model and if so, will be matched. - This error will be added to the pixel error around the position of point sources as follows: - sigma^2_i += 'psf_error_map'_j * (point_source_flux_i)**2 + Size can be larger or smaller than the pixel-sized PSF model and if so, will be matched. + This error will be added to the pixel error around the position of point sources as follows: + sigma^2_i += 'psf_error_map'_j * (point_source_flux_i)**2 :param point_source_supersampling_factor: int, supersampling factor of kernel_point_source. This is the input PSF to this class and does not need to be the choice in the modeling (thought preferred if modeling choses supersampling) diff --git a/lenstronomy/GalKin/analytic_kinematics.py b/lenstronomy/GalKin/analytic_kinematics.py index 4b8133e9f..798ba2be8 100644 --- a/lenstronomy/GalKin/analytic_kinematics.py +++ b/lenstronomy/GalKin/analytic_kinematics.py @@ -15,7 +15,7 @@ class AnalyticKinematics(Anisotropy): """ class to compute eqn 20 in Suyu+2010 with a Monte-Carlo from rendering from the - light profile distribution and displacing them with a Gaussian seeing convolution + light profile distribution and displacing them with a Gaussian seeing convolution. This class assumes spherical symmetry in light and mass distribution and - a Hernquist light profile (parameterised by the half-light radius) @@ -25,9 +25,7 @@ class to compute eqn 20 in Suyu+2010 with a Monte-Carlo from rendering from the the spectral rendering approach to compute the seeing convolved slit measurement is presented in Birrer et al. 2016. The stellar anisotropy is parameterised based on Osipkov 1979; Merritt 1985. - Units - ----- - all units are meant to be in angular arc seconds. The physical units are fold in through the angular diameter + All units are meant to be in angular arc seconds. The physical units are fold in through the angular diameter distances """ diff --git a/lenstronomy/GalKin/anisotropy.py b/lenstronomy/GalKin/anisotropy.py index 4d4e772dd..a2476e8c3 100644 --- a/lenstronomy/GalKin/anisotropy.py +++ b/lenstronomy/GalKin/anisotropy.py @@ -42,6 +42,7 @@ def __init__(self, anisotropy_type): def beta_r(self, r, **kwargs): """ returns the anisotropy parameter at a given radius + :param r: 3d radius :param kwargs: parameters of the specified anisotropy model :return: beta(r) @@ -73,7 +74,8 @@ def anisotropy_solution(self, r, **kwargs): def delete_anisotropy_cache(self): """ deletes cached interpolations for a fixed anisotropy model - :return: + + :return: None """ if hasattr(self._model, 'delete_cache'): self._model.delete_cache() diff --git a/lenstronomy/GalKin/aperture.py b/lenstronomy/GalKin/aperture.py index e3669d476..bc9287e44 100644 --- a/lenstronomy/GalKin/aperture.py +++ b/lenstronomy/GalKin/aperture.py @@ -26,7 +26,7 @@ def __init__(self, aperture_type, **kwargs_aperture): :param aperture_type: string :param kwargs_aperture: keyword arguments reflecting the aperture type chosen. - We refer to the specific class instances for documentation. + We refer to the specific class instances for documentation. """ if aperture_type == 'slit': self._aperture = Slit(**kwargs_aperture) diff --git a/lenstronomy/GalKin/aperture_types.py b/lenstronomy/GalKin/aperture_types.py index 12ee0a341..ee75e6412 100644 --- a/lenstronomy/GalKin/aperture_types.py +++ b/lenstronomy/GalKin/aperture_types.py @@ -39,6 +39,7 @@ def aperture_select(self, ra, dec): def num_segments(self): """ number of segments with separate measurements of the velocity dispersion + :return: int """ return 1 @@ -101,6 +102,7 @@ def aperture_select(self, ra, dec): def num_segments(self): """ number of segments with separate measurements of the velocity dispersion + :return: int """ return 1 @@ -161,6 +163,7 @@ def aperture_select(self, ra, dec): def num_segments(self): """ number of segments with separate measurements of the velocity dispersion + :return: int """ return 1 @@ -196,7 +199,7 @@ def __init__(self, r_bins, center_ra=0, center_dec=0): """ :param r_bins: array of radial bins to average the dispersion spectra in ascending order. - It starts with the inner-most edge to the outermost edge. + It starts with the inner-most edge to the outermost edge. :param center_ra: center of the sphere :param center_dec: center of the sphere """ @@ -228,7 +231,7 @@ def shell_ifu_select(ra, dec, r_bin, center_ra=0, center_dec=0): :param ra: angular coordinate of photon/ray :param dec: angular coordinate of photon/ray :param r_bin: array of radial bins to average the dispersion spectra in ascending order. - It starts with the inner-most edge to the outermost edge. + It starts with the inner-most edge to the outermost edge. :param center_ra: center of the sphere :param center_dec: center of the sphere :return: boolean, True if within the radial range, False otherwise diff --git a/lenstronomy/GalKin/cosmo.py b/lenstronomy/GalKin/cosmo.py index a2a36f5bb..ca6ccbf93 100644 --- a/lenstronomy/GalKin/cosmo.py +++ b/lenstronomy/GalKin/cosmo.py @@ -25,8 +25,9 @@ def __init__(self, d_d, d_s, d_ds): def arcsec2phys_lens(self, theta): """ converts are seconds to physical units on the deflector + :param theta: angle observed on the sky in units of arc seconds - :return: pyhsical distance of the angle in units of Mpc + :return: physical distance of the angle in units of Mpc """ return theta * const.arcsec * self.dd diff --git a/lenstronomy/GalKin/numeric_kinematics.py b/lenstronomy/GalKin/numeric_kinematics.py index 9174376b6..153eef943 100644 --- a/lenstronomy/GalKin/numeric_kinematics.py +++ b/lenstronomy/GalKin/numeric_kinematics.py @@ -179,7 +179,7 @@ def delete_cache(self): def _I_R_sigma2(self, R, kwargs_mass, kwargs_light, kwargs_anisotropy): """ - equation A15 in Mamon&Lokas 2005 as a logarithmic numerical integral (if option is chosen) + equation A15 in Mamon & Lokas 2005 as a logarithmic numerical integral (if option is chosen) :param R: 2d projected radius (in angular units) :param kwargs_mass: mass model parameters (following lenstronomy lens model conventions) diff --git a/lenstronomy/ImSim/Numerics/adaptive_numerics.py b/lenstronomy/ImSim/Numerics/adaptive_numerics.py index 35ade94dc..ce568ca21 100644 --- a/lenstronomy/ImSim/Numerics/adaptive_numerics.py +++ b/lenstronomy/ImSim/Numerics/adaptive_numerics.py @@ -27,7 +27,7 @@ def __init__(self, kernel_super, supersampling_factor, conv_supersample_pixels, :param supersampling_factor: factor of supersampling relative to pixel grid :param conv_supersample_pixels: bool array same size as data, pixels to be convolved and their light to be blurred :param supersampling_kernel_size: number of pixels (in units of the image pixels) that are convolved with the - supersampled kernel + supersampled kernel :param compute_pixels: bool array of size of image, these pixels (if True) will get blurred light from other pixels :param nopython: bool, numba jit setting to use python or compiled. :param cache: bool, numba jit setting to use cache diff --git a/lenstronomy/ImSim/Numerics/convolution.py b/lenstronomy/ImSim/Numerics/convolution.py index 766c83958..29a9c9987 100644 --- a/lenstronomy/ImSim/Numerics/convolution.py +++ b/lenstronomy/ImSim/Numerics/convolution.py @@ -152,7 +152,7 @@ def _static_pre_compute(self, image): # in1, s1, in2, s2 = in2, s2, in1, s1 # Speed up FFT by padding to optimal size for FFTPACK - fshape = [fftpack.helper.next_fast_len(int(d)) for d in shape] + fshape = [fftpack.next_fast_len(int(d)) for d in shape] fslice = tuple([slice(0, int(sz)) for sz in shape]) # Pre-1.9 NumPy FFT routines are not threadsafe. For older NumPys, make # sure we only call rfftn/irfftn from one thread at a time. @@ -191,7 +191,7 @@ def __init__(self, kernel_supersampled, supersampling_factor, supersampling_kern :param kernel_supersampled: kernel in supersampled pixels :param supersampling_factor: supersampling factor relative to the image pixel grid :param supersampling_kernel_size: number of pixels (in units of the image pixels) that are convolved with the - supersampled kernel + supersampled kernel """ # n_high = len(kernel_supersampled) self._supersampling_factor = supersampling_factor @@ -253,7 +253,7 @@ def __init__(self, sigma_list, fraction_list, pixel_scale, supersampling_factor= :param fraction_list: fraction of flux to be convoled with each Gaussian kernel :param pixel_scale: scale of pixel width (to convert sigmas into units of pixels) :param truncation: float. Truncate the filter at this many standard deviations. - Default is 4.0. + Default is 4.0. """ self._num_gaussians = len(sigma_list) self._sigmas_scaled = np.array(sigma_list) / pixel_scale @@ -276,10 +276,10 @@ def convolution2d(self, image): image_conv = None for i in range(self._num_gaussians): if image_conv is None: - image_conv = ndimage.filters.gaussian_filter(image, self._sigmas_scaled[i], mode='nearest', + image_conv = ndimage.gaussian_filter(image, self._sigmas_scaled[i], mode='nearest', truncate=self._truncation) * self._fraction_list[i] else: - image_conv += ndimage.filters.gaussian_filter(image, self._sigmas_scaled[i], mode='nearest', + image_conv += ndimage.gaussian_filter(image, self._sigmas_scaled[i], mode='nearest', truncate=self._truncation) * self._fraction_list[i] return image_conv diff --git a/lenstronomy/ImSim/Numerics/grid.py b/lenstronomy/ImSim/Numerics/grid.py index 2192ef825..8f0ffa04a 100644 --- a/lenstronomy/ImSim/Numerics/grid.py +++ b/lenstronomy/ImSim/Numerics/grid.py @@ -24,7 +24,7 @@ def __init__(self, nx, ny, transform_pix2angle, ra_at_xy_0, dec_at_xy_0, supersa :param supersampling_indexes: bool array of shape nx x ny, corresponding to pixels being super_sampled :param supersampling_factor: int, factor (per axis) of super-sampling :param flux_evaluate_indexes: bool array of shape nx x ny, corresponding to pixels being evaluated - (for both low and high res). Default is None, replaced by setting all pixels to being evaluated. + (for both low and high res). Default is None, replaced by setting all pixels to being evaluated. """ super(AdaptiveGrid, self).__init__(transform_pix2angle, ra_at_xy_0, dec_at_xy_0) self._nx = nx @@ -171,7 +171,7 @@ def __init__(self, nx, ny, transform_pix2angle, ra_at_xy_0, dec_at_xy_0, supersa :param dec_at_xy_0: dec coordinate at pixel (0,0) :param supersampling_factor: int, factor (per axis) of super-sampling :param flux_evaluate_indexes: bool array of shape nx x ny, corresponding to pixels being evaluated - (for both low and high res). Default is None, replaced by setting all pixels to being evaluated. + (for both low and high res). Default is None, replaced by setting all pixels to being evaluated. """ super(RegularGrid, self).__init__(transform_pix2angle, ra_at_xy_0, dec_at_xy_0) self._supersampling_factor = supersampling_factor diff --git a/lenstronomy/ImSim/Numerics/numerics.py b/lenstronomy/ImSim/Numerics/numerics.py index ceb89d759..09a12fdbe 100644 --- a/lenstronomy/ImSim/Numerics/numerics.py +++ b/lenstronomy/ImSim/Numerics/numerics.py @@ -24,19 +24,19 @@ def __init__(self, pixel_grid, psf, supersampling_factor=1, compute_mode='regula :param compute_mode: options are: 'regular', 'adaptive' :param supersampling_factor: int, factor of higher resolution sub-pixel sampling of surface brightness :param supersampling_convolution: bool, if True, performs (part of) the convolution on the super-sampled - grid/pixels + grid/pixels :param supersampling_kernel_size: int (odd number), size (in regular pixel units) of the super-sampled - convolution + convolution :param flux_evaluate_indexes: boolean 2d array of size of image (or None, then initiated as gird of True's). - Pixels indicated with True will be used to perform the surface brightness computation (and possible lensing - ray-shooting). Pixels marked as False will be assigned a flux value of zero (or ignored in the adaptive - convolution) + Pixels indicated with True will be used to perform the surface brightness computation (and possible lensing + ray-shooting). Pixels marked as False will be assigned a flux value of zero (or ignored in the adaptive + convolution) :param supersampled_indexes: 2d boolean array (only used in mode='adaptive') of pixels to be supersampled (in - surface brightness and if supersampling_convolution=True also in convolution). All other pixels not set to =True - will not be super-sampled. + surface brightness and if supersampling_convolution=True also in convolution). All other pixels not set to =True + will not be super-sampled. :param compute_indexes: 2d boolean array (only used in compute_mode='adaptive'), marks pixel that the response after - convolution is computed (all others =0). This can be set to likelihood_mask in the Likelihood module for - consistency. + convolution is computed (all others =0). This can be set to likelihood_mask in the Likelihood module for + consistency. :param point_source_supersampling_factor: super-sampling resolution of the point source placing :param convolution_kernel_size: int, odd number, size of convolution kernel. If None, takes size of point_source_kernel :param convolution_type: string, 'fft', 'grid', 'fft_static' mode of 2d convolution diff --git a/lenstronomy/ImSim/de_lens.py b/lenstronomy/ImSim/de_lens.py index 33f4fcbe1..4844d9a99 100644 --- a/lenstronomy/ImSim/de_lens.py +++ b/lenstronomy/ImSim/de_lens.py @@ -11,6 +11,7 @@ def get_param_WLS(A, C_D_inv, d, inv_bool=True): """ returns the parameter values given + :param A: response matrix Nd x Ns (Nd = # data points, Ns = # parameters) :param C_D_inv: inverse covariance matrix of the data, Nd x Nd, diagonal form :param d: data array, 1-d Nd @@ -44,6 +45,7 @@ def get_param_WLS(A, C_D_inv, d, inv_bool=True): def marginalisation_const(M_inv): """ get marginalisation constant 1/2 log(M_beta) for flat priors + :param M_inv: 2D covariance matrix :return: float """ diff --git a/lenstronomy/ImSim/differential_extinction.py b/lenstronomy/ImSim/differential_extinction.py index f76faf880..4111f2a98 100644 --- a/lenstronomy/ImSim/differential_extinction.py +++ b/lenstronomy/ImSim/differential_extinction.py @@ -14,7 +14,7 @@ def __init__(self, optical_depth_model=None, tau0_index=0): """ :param optical_depth_model: list of strings naming the profiles (same convention as LightModel module) - describing the optical depth of the extinction + describing the optical depth of the extinction """ if optical_depth_model is None: optical_depth_model = [] diff --git a/lenstronomy/ImSim/image2source_mapping.py b/lenstronomy/ImSim/image2source_mapping.py index bef75906c..5cfc137ba 100644 --- a/lenstronomy/ImSim/image2source_mapping.py +++ b/lenstronomy/ImSim/image2source_mapping.py @@ -24,12 +24,13 @@ class Image2SourceMapping(object): def __init__(self, lensModel, sourceModel): """ - :param lensModel: lenstronomy LensModel() class instance - :param sourceModel: LightModel () class instance - The lightModel includes: - - source_scale_factor_list: list of floats corresponding to the rescaled deflection angles to the specific source - components. None indicates that the list will be set to 1, meaning a single source plane model (in single lens plane mode). - - source_redshift_list: list of redshifts of the light components (in multi lens plane mode) + :param lensModel: LensModel() class instance + :param sourceModel: LightModel() class instance. + + The lightModel includes: + + - source_scale_factor_list: list of floats corresponding to the rescaled deflection angles to the specific source components. None indicates that the list will be set to 1, meaning a single source plane model (in single lens plane mode). + - source_redshift_list: list of redshifts of the light components (in multi lens plane mode) """ self._lightModel = sourceModel self._lensModel = lensModel diff --git a/lenstronomy/LensModel/LightConeSim/light_cone.py b/lenstronomy/LensModel/LightConeSim/light_cone.py index 51a5e547d..97a2f6b31 100644 --- a/lenstronomy/LensModel/LightConeSim/light_cone.py +++ b/lenstronomy/LensModel/LightConeSim/light_cone.py @@ -25,8 +25,10 @@ class to perform multi-plane ray-tracing from convergence maps at different reds def __init__(self, mass_map_list, grid_spacing_list, redshift_list): """ - :param mass_map_list: 2d numpy array of mass map (in units Msol) + :param mass_map_list: 2d numpy array of mass map + (in units physical Solar masses enclosed in each pixel/gird point of the map) :param grid_spacing_list: list of grid spacing of the individual mass maps + in units of physical Mpc :param redshift_list: list of redshifts of the mass maps """ diff --git a/lenstronomy/LensModel/MultiPlane/multi_plane.py b/lenstronomy/LensModel/MultiPlane/multi_plane.py index 644939357..5f7486133 100644 --- a/lenstronomy/LensModel/MultiPlane/multi_plane.py +++ b/lenstronomy/LensModel/MultiPlane/multi_plane.py @@ -25,16 +25,16 @@ def __init__(self, z_source, lens_model_list, lens_redshift_list, cosmo=None, nu :param lens_redshift_list: list of floats with redshifts of the lens models indicated in lens_model_list :param cosmo: instance of astropy.cosmology :param numerical_alpha_class: an instance of a custom class for use in NumericalAlpha() lens model - (see documentation in Profiles/numerical_alpha) + (see documentation in Profiles/numerical_alpha) :param kwargs_interp: interpolation keyword arguments specifying the numerics. See description in the Interpolate() class. Only applicable for 'INTERPOL' and 'INTERPOL_SCALED' models. :param observed_convention_index: a list of indices, corresponding to the lens_model_list element with same - index, where the 'center_x' and 'center_y' kwargs correspond to observed (lensed) positions, not physical - positions. The code will compute the physical locations when performing computations + index, where the 'center_x' and 'center_y' kwargs correspond to observed (lensed) positions, not physical + positions. The code will compute the physical locations when performing computations :param ignore_observed_positions: bool, if True, will ignore the conversion between observed to physical - position of deflectors + position of deflectors :param z_source_convention: float, redshift of a source to define the reduced deflection angles of the lens - models. If None, 'z_source' is used. + models. If None, 'z_source' is used. """ if z_source_convention is None: diff --git a/lenstronomy/LensModel/Profiles/cored_steep_ellipsoid.py b/lenstronomy/LensModel/Profiles/cored_steep_ellipsoid.py index 65b59f970..7ecec92ef 100644 --- a/lenstronomy/LensModel/Profiles/cored_steep_ellipsoid.py +++ b/lenstronomy/LensModel/Profiles/cored_steep_ellipsoid.py @@ -5,12 +5,13 @@ from lenstronomy.Util import param_util from lenstronomy.Util import util -__all__ = ['CSE', 'CSEMajorAxis', 'CSEMajorAxisSet'] +__all__ = ['CSE', 'CSEMajorAxis', 'CSEMajorAxisSet', 'CSEProductAvg', 'CSEProductAvgSet'] class CSE(LensProfileBase): """ Cored steep ellipsoid (CSE) + :param axis: 'major' or 'product_avg' ; whether to evaluate corresponding to r= major axis or r= sqrt(ab) source: Keeton and Kochanek (1998) Oguri 2021: https://arxiv.org/pdf/2106.11464.pdf @@ -28,8 +29,13 @@ class CSE(LensProfileBase): lower_limit_default = {'A': -1000, 's': 0, 'e1': -0.5, 'e2': -0.5, 'center_x': -100, 'center_y': -100} upper_limit_default = {'A': 1000, 's': 10000, 'e1': 0.5, 'e2': 0.5, 'center_x': -100, 'center_y': -100} - def __init__(self): - self.major_axis_model = CSEMajorAxis() + def __init__(self, axis='product_avg'): + if axis == 'major': + self.major_axis_model = CSEMajorAxis() + elif axis == 'product_avg': + self.major_axis_model = CSEProductAvg() + else: + raise ValueError("axis must be set to 'major' or 'product_avg'. Input is %s ." % axis) super(CSE, self).__init__() def function(self, x, y, a, s, e1, e2, center_x, center_y): @@ -251,3 +257,143 @@ def hessian(self, x, y, a_list, s_list, q): f_xy += f_xy_ f_yy += f_yy_ return f_xx, f_xy, f_xy, f_yy + + +class CSEProductAvg(LensProfileBase): + """ + Cored steep ellipsoid (CSE) evaluated at the product-averaged radius sqrt(ab), + such that mass is not changed when increasing ellipticity + + Same as CSEMajorAxis but evaluated at r=sqrt(q)*r_original + + Keeton and Kochanek (1998) + Oguri 2021: https://arxiv.org/pdf/2106.11464.pdf + + .. math:: + \\kappa(u;s) = \\frac{A}{2(s^2 + \\xi^2)^{3/2}} + + with + + .. math:: + \\xi(x, y) = \\sqrt{qx^2 + \\frac{y^2}{q}} + + + + """ + param_names = ['A', 's', 'q', 'center_x', 'center_y'] + lower_limit_default = {'A': -1000, 's': 0, 'q': 0.001, 'center_x': -100, 'center_y': -100} + upper_limit_default = {'A': 1000, 's': 10000, 'q': 0.99999, 'e2': 0.5, 'center_x': -100, 'center_y': -100} + + def __init__(self): + super(CSEProductAvg, self).__init__() + self.MA_class = CSEMajorAxis() + + @staticmethod + def _convert2prodavg(x, y, a, s, q): + """ + converts coordinates and re-normalizes major-axis parameterization to instead be wrt. product-averaged + """ + a = a / q + x = x * np.sqrt(q) + y = y * np.sqrt(q) + return x, y, a, s, q + + def function(self, x, y, a, s, q): + """ + :param x: coordinate in image plane (angle) + :param y: coordinate in image plane (angle) + :param a: lensing strength + :param s: core radius + :param q: axis ratio + :return: lensing potential + """ + x, y, a, s, q = self._convert2prodavg(x, y, a, s, q) + return self.MA_class.function(x, y, a, s, q) + + def derivatives(self, x, y, a, s, q): + """ + :param x: coordinate in image plane (angle) + :param y: coordinate in image plane (angle) + :param a: lensing strength + :param s: core radius + :param q: axis ratio + :return: deflection in x- and y-direction + """ + x, y, a, s, q = self._convert2prodavg(x, y, a, s, q) + af_x, af_y = self.MA_class.derivatives(x, y, a, s, q) + # extra sqrt(q) factor from taking derivative of transformed coordinate + return np.sqrt(q) * af_x, np.sqrt(q) * af_y + + def hessian(self, x, y, a, s, q): + """ + :param x: coordinate in image plane (angle) + :param y: coordinate in image plane (angle) + :param a: lensing strength + :param s: core radius + :param q: axis ratio + :return: hessian elements f_xx, f_xy, f_yx, f_yy + """ + x, y, a, s, q = self._convert2prodavg(x, y, a, s, q) + af_xx, af_xy, af_xy, af_yy = self.MA_class.hessian(x, y, a, s, q) + # two sqrt(q) factors from taking derivatives of transformed coordinate + return q * af_xx, q * af_xy, q * af_xy, q * af_yy + + +class CSEProductAvgSet(LensProfileBase): + """ + a set of CSE profiles along a joint center and axis + """ + + def __init__(self): + self.major_axis_model = CSEProductAvg() + super(CSEProductAvgSet, self).__init__() + + def function(self, x, y, a_list, s_list, q): + """ + + :param x: coordinate in image plane (angle) + :param y: coordinate in image plane (angle) + :param a_list: list of lensing strength + :param s_list: list of core radius + :param q: axis ratio + :return: lensing potential + """ + f_ = np.zeros_like(x) + for a, s in zip(a_list, s_list): + f_ += self.major_axis_model.function(x, y, a, s, q) + return f_ + + def derivatives(self, x, y, a_list, s_list, q): + """ + + :param x: coordinate in image plane (angle) + :param y: coordinate in image plane (angle) + :param a_list: list of lensing strength + :param s_list: list of core radius + :param q: axis ratio + :return: deflection in x- and y-direction + """ + f_x, f_y = np.zeros_like(x), np.zeros_like(y) + for a, s in zip(a_list, s_list): + f_x_, f_y_ = self.major_axis_model.derivatives(x, y, a, s, q) + f_x += f_x_ + f_y += f_y_ + return f_x, f_y + + def hessian(self, x, y, a_list, s_list, q): + """ + + :param x: coordinate in image plane (angle) + :param y: coordinate in image plane (angle) + :param a_list: list of lensing strength + :param s_list: list of core radius + :param q: axis ratio + :return: hessian elements f_xx, f_xy, f_yx, f_yy + """ + f_xx, f_xy, f_yy = np.zeros_like(x), np.zeros_like(x), np.zeros_like(x) + for a, s in zip(a_list, s_list): + f_xx_, f_xy_, _, f_yy_ = self.major_axis_model.hessian(x, y, a, s, q) + f_xx += f_xx_ + f_xy += f_xy_ + f_yy += f_yy_ + return f_xx, f_xy, f_xy, f_yy diff --git a/lenstronomy/LensModel/Profiles/curved_arc_tan_diff.py b/lenstronomy/LensModel/Profiles/curved_arc_tan_diff.py index be06a413f..626abf2ac 100644 --- a/lenstronomy/LensModel/Profiles/curved_arc_tan_diff.py +++ b/lenstronomy/LensModel/Profiles/curved_arc_tan_diff.py @@ -26,7 +26,7 @@ class CurvedArcTanDiff(LensProfileBase): """ param_names = ['tangential_stretch', 'radial_stretch', 'curvature', 'dtan_dtan', 'direction', 'center_x', 'center_y'] lower_limit_default = {'tangential_stretch': -100, 'radial_stretch': -5, 'curvature': 0.000001, 'dtan_dtan': -10, 'direction': -np.pi, 'center_x': -100, 'center_y': -100} - upper_limit_default = {'tangential_stretch': 100, 'radial_stretch': 5, 'curvature': 100, 'dtan_dtab': 10, 'direction': np.pi, 'center_x': 100, 'center_y': 100} + upper_limit_default = {'tangential_stretch': 100, 'radial_stretch': 5, 'curvature': 100, 'dtan_dtan': 10, 'direction': np.pi, 'center_x': 100, 'center_y': 100} def __init__(self): self._sie = SIE(NIE=True) diff --git a/lenstronomy/LensModel/Profiles/epl.py b/lenstronomy/LensModel/Profiles/epl.py index 15499ab7b..f33063b2f 100644 --- a/lenstronomy/LensModel/Profiles/epl.py +++ b/lenstronomy/LensModel/Profiles/epl.py @@ -26,7 +26,7 @@ class EPL(LensProfileBase): In terms of eccentricities, this profile is defined as .. math:: - \\kappa(r) = \\frac{3-\\gamma}{2} \\left(\\frac{\\theta'_{E}}{r \\sqrt{1 − e*\\cos(2*\\phi)}} \\right)^{\\gamma-1} + \\kappa(r) = \\frac{3-\\gamma}{2} \\left(\\frac{\\theta'_{E}}{r \\sqrt{1 - e*\\cos(2*\\phi)}} \\right)^{\\gamma-1} with :math:`\\epsilon` is the ellipticity defined as diff --git a/lenstronomy/LensModel/Profiles/general_nfw.py b/lenstronomy/LensModel/Profiles/general_nfw.py index ffad78dab..5006d58d8 100644 --- a/lenstronomy/LensModel/Profiles/general_nfw.py +++ b/lenstronomy/LensModel/Profiles/general_nfw.py @@ -14,13 +14,11 @@ class GNFW(LensProfileBase): This class contains a double power law profile with flexible inner and outer logarithmic slopes g and n .. math:: - \\rho = \\rho_0 (x^g (1+x^2)^((n-g)/2))^{-1} + \\rho(r) = \\frac{\\rho_0}{r^{\\gamma}} \\frac{Rs^{n}}{\\left(r^2 + Rs^2 \\right)^{(n - \\gamma)/2}} For g = 1.0 and n=3, it is approximately the same as an NFW profile - The original reference is [1]_ + The original reference is [1]_. - References - ---------- .. [1] Munoz, Kochanek and Keeton, (2001), astro-ph/0103009, doi:10.1086/322314 TODO: implement the gravitational potential for this profile @@ -70,14 +68,15 @@ def hessian(self, x, y, Rs, alpha_Rs, gamma_inner, gamma_outer, center_x=0, cent y_ = y - center_y R = np.sqrt(x_ ** 2 + y_ ** 2) R = np.maximum(R, 0.00000001) - kappa = self.density_2d(R, 0, Rs, rho0_input,gamma_inner, gamma_outer) - gamma1, gamma2 = self.nfwGamma(R, Rs, rho0_input,gamma_inner, gamma_outer, x_, y_) + kappa = self.density_2d(R, 0, Rs, rho0_input, gamma_inner, gamma_outer) + gamma1, gamma2 = self.nfwGamma(R, Rs, rho0_input, gamma_inner, gamma_outer, x_, y_) f_xx = kappa + gamma1 f_yy = kappa - gamma1 f_xy = gamma2 return f_xx, f_xy, f_xy, f_yy - def density(self, R, Rs, rho0, gamma_inner, gamma_outer): + @staticmethod + def density(R, Rs, rho0, gamma_inner, gamma_outer): """ three dimensional NFW profile @@ -90,7 +89,7 @@ def density(self, R, Rs, rho0, gamma_inner, gamma_outer): """ x = R/Rs outer_slope = (gamma_outer-gamma_inner)/2 - return rho0 / (x**gamma_inner * (1 +x ** 2) ** outer_slope) + return rho0 / (x**gamma_inner * (1 + x ** 2) ** outer_slope) def density_lens(self, r, Rs, alpha_Rs, gamma_inner, gamma_outer): """ @@ -114,7 +113,7 @@ def density_2d(self, x, y, Rs, rho0, gamma_inner, gamma_outer, center_x=0, cente :param x: angular position (normally in units of arc seconds) :param y: angular position (normally in units of arc seconds) :param Rs: turn over point in the slope of the NFW profile in angular unit - :param alpha_Rs: deflection (angular units) at projected Rs + :param rho0: density normalization at Rs :param gamma_inner: logarithmic profile slope interior to Rs :param gamma_outer: logarithmic profile slope outside Rs :param center_x: profile center (same units as x) @@ -128,7 +127,8 @@ def density_2d(self, x, y, Rs, rho0, gamma_inner, gamma_outer, center_x=0, cente Fx = self._f(x, gamma_inner, gamma_outer) return 2 * rho0 * Rs * Fx - def mass_3d(self, r, Rs, rho0, gamma_inner, gamma_outer): + @staticmethod + def mass_3d(r, Rs, rho0, gamma_inner, gamma_outer): """ mass enclosed a 3d sphere or radius r @@ -166,7 +166,7 @@ def mass_2d(self, R, Rs, rho0, gamma_inner, gamma_outer): """ mass enclosed a 2d cylinder or projected radius R - :param r: 3d radius + :param R: 3d radius :param Rs: scale radius :param rho0: central density normalization :param gamma_inner: logarithmic profile slope interior to Rs @@ -184,7 +184,7 @@ def nfwAlpha(self, R, Rs, rho0, gamma_inner, gamma_outer, ax_x, ax_y): deflection angel of NFW profile (times Sigma_crit D_OL) along the projection to coordinate 'axis' - :param r: 3d radius + :param R: 3d radius :param Rs: scale radius :param rho0: central density normalization :param gamma_inner: logarithmic profile slope interior to Rs @@ -204,7 +204,7 @@ def nfwGamma(self, R, Rs, rho0, gamma_inner, gamma_outer, ax_x, ax_y): shear gamma of NFW profile (times Sigma_crit) along the projection to coordinate 'axis' - :param r: 3d radius + :param R: 3d radius :param Rs: scale radius :param rho0: central density normalization :param gamma_inner: logarithmic profile slope interior to Rs @@ -231,8 +231,8 @@ def _f(X, g, n): :param n: logarithmic profile slope exterior to Rs :return: solution to the projection integral """ - if n==3: - n=3.001 # for numerical stability + if n == 3: + n = 3.001 # for numerical stability hyp2f1_term = hyp2f1((n-1)/2, g/2, n/2, 1/(1+X**2)) beta_term = beta((n-1)/2, 0.5) return 0.5 * beta_term * hyp2f1_term * (1+X**2) ** ((1-n)/2) @@ -248,8 +248,8 @@ def _g(X, g, n): :param n: logarithmic profile slope exterior to Rs :return: solution of the integral over projected mass """ - if n==3: - n=3.001 # for numerical stability + if n == 3: + n = 3.001 # for numerical stability xi = 1 + X**2 hyp2f1_term = hyp2f1((n - 3) / 2, g / 2, n / 2, 1 / xi) beta_term_1 = beta((n - 3) / 2, (3-g)/2) @@ -269,7 +269,7 @@ def alpha2rho0(self, alpha_Rs, Rs, gamma_inner, gamma_outer): """ gx = self._g(1.0, gamma_inner, gamma_outer) - rho0 = alpha_Rs / (4. * Rs **2 * gx / 1.0 ** 2) + rho0 = alpha_Rs / (4. * Rs ** 2 * gx / 1.0 ** 2) return rho0 def rho02alpha(self, rho0, Rs, gamma_inner, gamma_outer): @@ -284,5 +284,5 @@ def rho02alpha(self, rho0, Rs, gamma_inner, gamma_outer): :return: deflection angle at RS """ gx = self._g(1.0, gamma_inner, gamma_outer) - alpha_Rs = rho0 * (4. * Rs **2 * gx / 1.0 ** 2) + alpha_Rs = rho0 * (4. * Rs ** 2 * gx / 1.0 ** 2) return alpha_Rs diff --git a/lenstronomy/LensModel/Profiles/nfw.py b/lenstronomy/LensModel/Profiles/nfw.py index 8237ad422..7d57ba70d 100644 --- a/lenstronomy/LensModel/Profiles/nfw.py +++ b/lenstronomy/LensModel/Profiles/nfw.py @@ -19,6 +19,7 @@ class NFW(LensProfileBase): Examples for converting angular to physical mass units ------------------------------------------------------ + >>> from lenstronomy.Cosmo.lens_cosmo import LensCosmo >>> from astropy.cosmology import FlatLambdaCDM >>> cosmo = FlatLambdaCDM(H0=70, Om0=0.3, Ob0=0.05) @@ -34,7 +35,7 @@ class NFW(LensProfileBase): The lens model calculation uses angular units as arguments! So to execute a deflection angle calculation one uses - >>> from lenstronomy.LensModel.Profiles.hernquist import NFW + >>> from lenstronomy.LensModel.Profiles.nfw import NFW >>> nfw = NFW() >>> alpha_x, alpha_y = nfw.derivatives(x=1, y=1, Rs=Rs_angle, alpha_Rs=alpha_Rs, center_x=0, center_y=0) diff --git a/lenstronomy/LensModel/Profiles/nfw_ellipse_cse.py b/lenstronomy/LensModel/Profiles/nfw_ellipse_cse.py index 5928e8143..98e74375f 100644 --- a/lenstronomy/LensModel/Profiles/nfw_ellipse_cse.py +++ b/lenstronomy/LensModel/Profiles/nfw_ellipse_cse.py @@ -4,7 +4,7 @@ from lenstronomy.Util import util from lenstronomy.LensModel.Profiles.nfw import NFW from lenstronomy.LensModel.Profiles.nfw_ellipse import NFW_ELLIPSE -from lenstronomy.LensModel.Profiles.cored_steep_ellipsoid import CSEMajorAxisSet +from lenstronomy.LensModel.Profiles.cored_steep_ellipsoid import CSEProductAvgSet import lenstronomy.Util.param_util as param_util __all__ = ['NFW_ELLIPSE_CSE'] @@ -15,6 +15,7 @@ class NFW_ELLIPSE_CSE(NFW_ELLIPSE): this class contains functions concerning the NFW profile with an ellipticity defined in the convergence parameterization of alpha_Rs and Rs is the same as for the spherical NFW profile Approximation with CSE profile introduced by Oguri 2021: https://arxiv.org/pdf/2106.11464.pdf + Match to NFW using CSEs is approximate: kappa matches to ~1-2% relation are: R_200 = c * Rs @@ -28,9 +29,10 @@ class NFW_ELLIPSE_CSE(NFW_ELLIPSE): def __init__(self, high_accuracy=True): """ - :param high_accuracy: boolean, if True uses a more accurate larger set of CSE profiles (see Oguri 2021) + :param high_accuracy: if True uses a more accurate larger set of CSE profiles (see Oguri 2021) + :type high_accuracy: boolean """ - self.cse_major_axis_set = CSEMajorAxisSet() + self.cse_major_axis_set = CSEProductAvgSet() self.nfw = NFW() if high_accuracy is True: # Table 1 in Oguri 2021 @@ -167,6 +169,6 @@ def _normalization(self, alpha_Rs, Rs, q): :return: normalization (m) """ rho0 = self.nfw.alpha2rho0(alpha_Rs, Rs) - rs_ = Rs / np.sqrt(q) + rs_ = Rs const = 4 * rho0 * rs_ ** 3 return const diff --git a/lenstronomy/LensModel/Profiles/sersic.py b/lenstronomy/LensModel/Profiles/sersic.py index 064e4cd22..ae8ed4a56 100644 --- a/lenstronomy/LensModel/Profiles/sersic.py +++ b/lenstronomy/LensModel/Profiles/sersic.py @@ -1,5 +1,4 @@ __author__ = 'sibirrer' -#this file contains a class to make a gaussian import numpy as np import lenstronomy.Util.util as util @@ -12,6 +11,43 @@ class Sersic(SersicUtil, LensProfileBase): """ this class contains functions to evaluate a Sersic mass profile: https://arxiv.org/pdf/astro-ph/0311559.pdf + + .. math:: + \\kappa(R) = \\kappa_{\\rm eff} \\exp \\left[ -b_n (R/R_{\\rm Sersic})^{\\frac{1}{n}}\\right] + + with :math:`b_{n}\\approx 1.999n-0.327` + + Examples + -------- + + Example for converting physical mass units into convergence units used in the definition of this profile. + + We first define an AstroPy cosmology instance and a LensCosmo class instance with a lens and source redshift. + + >>> from lenstronomy.Cosmo.lens_cosmo import LensCosmo + >>> from astropy.cosmology import FlatLambdaCDM + >>> cosmo = FlatLambdaCDM(H0=70, Om0=0.3, Ob0=0.05) + >>> lens_cosmo = LensCosmo(z_lens=0.5, z_source=1.5, cosmo=cosmo) + + We define the half-light radius R_sersic (arc seconds on the sky) and Sersic index n_sersic + + >>> R_sersic = 2 + >>> n_sersic = 4 + + Here we compute k_eff, the convergence at the half-light radius R_sersic for a stellar mass in Msun + + >>> k_eff = lens_cosmo.sersic_m_star2k_eff(m_star=10**11.5, R_sersic=R_sersic, n_sersic=n_sersic) + + And here we perform the inverse calculation given k_eff to return the physical stellar mass. + + >>> m_star = lens_cosmo.sersic_k_eff2m_star(k_eff=k_eff, R_sersic=R_sersic, n_sersic=n_sersic) + + The lens model calculation uses angular units as arguments! So to execute a deflection angle calculation one uses + + >>> from lenstronomy.LensModel.Profiles.sersic import Sersic + >>> sersic = Sersic() + >>> alpha_x, alpha_y = sersic.derivatives(x=1, y=1, k_eff=k_eff, R_sersic=R_sersic, center_x=0, center_y=0) + """ param_names = ['k_eff', 'R_sersic', 'n_sersic', 'center_x', 'center_y'] lower_limit_default = {'k_eff': 0, 'R_sersic': 0, 'n_sersic': 0.5, 'center_x': -100, 'center_y': -100} @@ -69,10 +105,6 @@ def hessian(self, x, y, n_sersic, R_sersic, k_eff, center_x=0, center_y=0): d_alpha_dr = self.d_alpha_dr(x, y, n_sersic, R_sersic, k_eff, center_x, center_y) alpha = -self.alpha_abs(x, y, n_sersic, R_sersic, k_eff, center_x, center_y) - #f_xx_ = d_alpha_dr * calc_util.d_r_dx(x_, y_) * x_/r + alpha * calc_util.d_x_diffr_dx(x_, y_) - #f_yy_ = d_alpha_dr * calc_util.d_r_dy(x_, y_) * y_/r + alpha * calc_util.d_y_diffr_dy(x_, y_) - #f_xy_ = d_alpha_dr * calc_util.d_r_dy(x_, y_) * x_/r + alpha * calc_util.d_x_diffr_dy(x_, y_) - f_xx = -(d_alpha_dr/r + alpha/r**2) * x_**2/r + alpha/r f_yy = -(d_alpha_dr/r + alpha/r**2) * y_**2/r + alpha/r f_xy = -(d_alpha_dr/r + alpha/r**2) * x_*y_/r diff --git a/lenstronomy/LensModel/Profiles/shear.py b/lenstronomy/LensModel/Profiles/shear.py index 031fb1ba9..75a607fc4 100644 --- a/lenstronomy/LensModel/Profiles/shear.py +++ b/lenstronomy/LensModel/Profiles/shear.py @@ -75,8 +75,8 @@ class to model a shear field with shear strength and direction. The translation is as follow: .. math:: - \\gamma_1 = \\gamma_{ext} \\cos(2 \\phi_{ext} - \\gamma_2 = \\gamma_{ext} \\sin(2 \\phi_{ext} + \\gamma_1 = \\gamma_{ext} \\cos(2 \\phi_{ext}) + \\gamma_2 = \\gamma_{ext} \\sin(2 \\phi_{ext}) """ param_names = ['gamma_ext', 'psi_ext', 'ra_0', 'dec_0'] @@ -121,7 +121,7 @@ class ShearReduced(LensProfileBase): To keep the magnification at unity, it requires .. math:: - (1-\\kappa)^2 - \\gamma_1^2 - \\gamma_2^ = 1 + (1-\\kappa)^2) - \\gamma_1^2 - \\gamma_2^ = 1 Thus, for given pair of reduced shear :math:`(\\gamma'_1, \\gamma'_2)`, an additional convergence term is calculated and added to the lensing distortions. diff --git a/lenstronomy/LensModel/Profiles/sie.py b/lenstronomy/LensModel/Profiles/sie.py index d2dcf66d3..ef2cc3932 100644 --- a/lenstronomy/LensModel/Profiles/sie.py +++ b/lenstronomy/LensModel/Profiles/sie.py @@ -106,6 +106,7 @@ def hessian(self, x, y, theta_E, e1, e2, center_x=0, center_y=0): def theta2rho(theta_E): """ converts projected density parameter (in units of deflection) into 3d density parameter + :param theta_E: :return: """ @@ -117,6 +118,7 @@ def theta2rho(theta_E): def mass_3d(r, rho0, e1=0, e2=0): """ mass enclosed a 3d sphere or radius r + :param r: radius in angular units :param rho0: density at angle=1 :return: mass in angular units @@ -138,6 +140,7 @@ def mass_3d_lens(self, r, theta_E, e1=0, e2=0): def mass_2d(self, r, rho0, e1=0, e2=0): """ mass enclosed projected 2d sphere of radius r + :param r: :param rho0: :param e1: @@ -163,6 +166,7 @@ def mass_2d_lens(self, r, theta_E, e1=0, e2=0): def grav_pot(self, x, y, rho0, e1=0, e2=0, center_x=0, center_y=0): """ gravitational potential (modulo 4 pi G and rho0 in appropriate units) + :param x: :param y: :param rho0: @@ -197,6 +201,7 @@ def density_lens(self, r, theta_E, e1=0, e2=0): def density(r, rho0, e1=0, e2=0): """ computes the density + :param r: radius in angles :param rho0: density at angle=1 :return: density at r @@ -208,6 +213,7 @@ def density(r, rho0, e1=0, e2=0): def density_2d(x, y, rho0, e1=0, e2=0, center_x=0, center_y=0): """ projected density + :param x: :param y: :param rho0: diff --git a/lenstronomy/LensModel/Profiles/spemd.py b/lenstronomy/LensModel/Profiles/spemd.py index 088920be6..1dfe8deb3 100644 --- a/lenstronomy/LensModel/Profiles/spemd.py +++ b/lenstronomy/LensModel/Profiles/spemd.py @@ -216,6 +216,7 @@ def convert_params(theta_E, gamma, q, s_scale): def is_not_empty(x1, x2): """ Check if float or not an empty array + :return: True if x1 and x2 are either floats/ints or an non-empty array, False if e.g. objects are [] :rtype: bool """ diff --git a/lenstronomy/LensModel/Profiles/spep.py b/lenstronomy/LensModel/Profiles/spep.py index b4847a208..4d579cd7a 100644 --- a/lenstronomy/LensModel/Profiles/spep.py +++ b/lenstronomy/LensModel/Profiles/spep.py @@ -30,10 +30,10 @@ def function(self, x, y, theta_E, gamma, e1, e2, center_x=0, center_y=0): :type theta_E: float. :param gamma: power law slope of mass profifle :type gamma: <2 float - :param q: Axis ratio - :type q: 0 infinity :return: mass inside radius r @@ -178,6 +184,7 @@ def mass_3d(self, r, rho0, r_core, gamma): def mass_3d_lens(self, r, sigma0, r_core, gamma): """ mass enclosed a 3d sphere or radius r + :param r: radius [arcsec] :param sigma0: convergence at r = 0 :param r_core: core radius [arcsec] @@ -190,9 +197,10 @@ def mass_3d_lens(self, r, sigma0, r_core, gamma): def mass_2d(self, r, rho0, r_core, gamma): """ mass enclosed projected 2d disk of radius r + :param r: radius [arcsec] :param rho0: density at r = 0 in units [rho_0_physical / sigma_crit] (which should be equal to [1/arcsec]) - where rho_0_physical is a physical density normalization and sigma_crit is the critical density for lensing + where rho_0_physical is a physical density normalization and sigma_crit is the critical density for lensing :param r_core: core radius [arcsec] :param gamma: logarithmic slope at r -> infinity :return: projected mass inside disk of radius r @@ -202,9 +210,10 @@ def mass_2d(self, r, rho0, r_core, gamma): def mass_2d_lens(self, r, sigma0, r_core, gamma): """ mass enclosed projected 2d disk of radius r + :param r: radius [arcsec] :param sigma0: convergence at r = 0 - where rho_0_physical is a physical density normalization and sigma_crit is the critical density for lensing + where rho_0_physical is a physical density normalization and sigma_crit is the critical density for lensing :param r_core: core radius [arcsec] :param gamma: logarithmic slope at r -> infinity :return: projected mass inside disk of radius r @@ -216,6 +225,7 @@ def mass_2d_lens(self, r, sigma0, r_core, gamma): def _safe_r_division(r, r_core, x_min=1e-6): """ Avoids accidental division by 0 + :param r: radius in arcsec :param r_core: core radius in arcsec :return: a minimum value of r @@ -231,6 +241,7 @@ def _safe_r_division(r, r_core, x_min=1e-6): def _sigma2rho0(sigma0, r_core): """ Converts the convergence normalization to the 3d normalization + :param sigma0: convergence at r=0 :param r_core: core radius [arcsec] :return: density normalization in units 1/arcsec, or rho_0_physical / sigma_crit @@ -241,6 +252,7 @@ def _sigma2rho0(sigma0, r_core): def _rho02sigma(rho0, r_core): """ Converts the convergence normalization to the 3d normalization + :param rho0: convergence at r=0 :param r_core: core radius [arcsec] :return: density normalization in units 1/arcsec, or rho_0_physical / sigma_crit diff --git a/lenstronomy/LensModel/Profiles/spp.py b/lenstronomy/LensModel/Profiles/spp.py index 0834e604e..9927ea929 100644 --- a/lenstronomy/LensModel/Profiles/spp.py +++ b/lenstronomy/LensModel/Profiles/spp.py @@ -20,6 +20,8 @@ def function(self, x, y, theta_E, gamma, center_x=0, center_y=0): """ :param x: set of x-coordinates :type x: array of size (n) + :param y: set of y-coordinates + :type y: array of size (n) :param theta_E: Einstein radius of lens :type theta_E: float. :param gamma: power law slope of mass profile @@ -85,6 +87,7 @@ def hessian(self, x, y, theta_E, gamma, center_x=0., center_y=0.): def rho2theta(rho0, gamma): """ converts 3d density into 2d projected density parameter + :param rho0: :param gamma: :return: @@ -99,6 +102,7 @@ def rho2theta(rho0, gamma): def theta2rho(theta_E, gamma): """ converts projected density parameter (in units of deflection) into 3d density parameter + :param theta_E: :param gamma: :return: @@ -112,9 +116,10 @@ def theta2rho(theta_E, gamma): def mass_3d(r, rho0, gamma): """ mass enclosed a 3d sphere or radius r + :param r: - :param a: - :param s: + :param rho0: + :param gamma: :return: """ mass_3d = 4 * np.pi * rho0 /(-gamma + 3) * r ** (-gamma + 3) @@ -134,10 +139,10 @@ def mass_3d_lens(self, r, theta_E, gamma): def mass_2d(self, r, rho0, gamma): """ mass enclosed projected 2d sphere of radius r + :param r: :param rho0: - :param a: - :param s: + :param gamma: :return: """ alpha = np.sqrt(np.pi) * special.gamma(1. / 2 * (-1 + gamma)) / special.gamma(gamma / 2.) * r ** (2 - gamma)/(3 - gamma) * 2 * rho0 @@ -158,11 +163,11 @@ def mass_2d_lens(self, r, theta_E, gamma): def grav_pot(self, x, y, rho0, gamma, center_x=0, center_y=0): """ gravitational potential (modulo 4 pi G and rho0 in appropriate units) + :param x: :param y: :param rho0: - :param a: - :param s: + :param gamma: :param center_x: :param center_y: :return: @@ -178,11 +183,10 @@ def grav_pot(self, x, y, rho0, gamma, center_x=0, center_y=0): def density(r, rho0, gamma): """ computes the density - :param x: - :param y: + + :param r: :param rho0: - :param a: - :param s: + :param gamma: :return: """ rho = rho0 / r**gamma @@ -201,11 +205,11 @@ def density_lens(self, r, theta_E, gamma): def density_2d(x, y, rho0, gamma, center_x=0, center_y=0): """ projected density + :param x: :param y: :param rho0: - :param a: - :param s: + :param gamma: :param center_x: :param center_y: :return: diff --git a/lenstronomy/LensModel/Profiles/tnfw.py b/lenstronomy/LensModel/Profiles/tnfw.py index f252c2388..3aea301ca 100644 --- a/lenstronomy/LensModel/Profiles/tnfw.py +++ b/lenstronomy/LensModel/Profiles/tnfw.py @@ -58,6 +58,7 @@ def function(self, x, y, Rs, alpha_Rs, r_trunc, center_x=0, center_y=0): def _L(self, x, tau): """ Logarithm that appears frequently + :param x: r/Rs :param tau: t/Rs :return: @@ -68,6 +69,7 @@ def _L(self, x, tau): def F(self, x): """ Classic NFW function in terms of arctanh and arctan + :param x: r/Rs :return: """ diff --git a/lenstronomy/LensModel/Profiles/uldm.py b/lenstronomy/LensModel/Profiles/uldm.py index 44ff9ded4..f47d8d1eb 100644 --- a/lenstronomy/LensModel/Profiles/uldm.py +++ b/lenstronomy/LensModel/Profiles/uldm.py @@ -2,10 +2,9 @@ # this file contains a class to compute the Ultra Light Dark Matter soliton profile import numpy as np -import scipy.interpolate as interp from scipy.special import gamma, hyp2f1 from lenstronomy.LensModel.Profiles.base_profile import LensProfileBase -import lenstronomy.Util.constants as const + __all__ = ['Uldm'] @@ -29,24 +28,27 @@ class Uldm(LensProfileBase): different from 8 to model solitons which feel the influence of background potential (see 2105.10873) The profile has, as parameters: - :param kappa_0: central convergence - :param theta_c: core radius (in arcseconds) - :param slope: exponent entering the profile, default value is 8 + + - kappa_0: central convergence + - theta_c: core radius (in arcseconds) + - slope: exponent entering the profile, default value is 8 """ _s = 0.000001 # numerical limit for minimal radius param_names = ['kappa_0', 'theta_c', 'slope', 'center_x', 'center_y'] lower_limit_default = {'kappa_0': 0, 'theta_c': 0, 'slope': 3.5, 'center_x': -100, 'center_y': -100} upper_limit_default = {'kappa_0': 1., 'theta_c': 100, 'slope': 10, 'center_x': 100, 'center_y': 100} - def rhotilde(self, kappa_0, theta_c, slope=8): + @staticmethod + def rhotilde(kappa_0, theta_c, slope=8): """ Computes the central density in angular units + :param kappa_0: central convergence of profile :param theta_c: core radius (in arcsec) :param slope: exponent entering the profile :return: central density in 1/arcsec """ - a_factor_sqrt = np.sqrt( (0.5)**(-1/slope) -1) + a_factor_sqrt = np.sqrt(0.5**(-1/slope) - 1) num_factor = gamma(slope) / gamma(slope - 1/2) * a_factor_sqrt / np.sqrt(np.pi) return kappa_0 * num_factor / theta_c @@ -76,7 +78,8 @@ def function(self, x, y, kappa_0, theta_c, center_x=0, center_y=0, slope=8): hyp3f2(1, 1, slope - 0.5, 2, 2, -(a_factor_sqrt * r_i / theta_c)**2.) for r_i in r], dtype=float) return hypgeom - def alpha_radial(self, r, kappa_0, theta_c, slope=8): + @staticmethod + def alpha_radial(r, kappa_0, theta_c, slope=8): """ returns the radial part of the deflection angle @@ -86,8 +89,8 @@ def alpha_radial(self, r, kappa_0, theta_c, slope=8): :param r: radius where the deflection angle is computed :return: radial deflection angle """ - a_factor = (0.5)**(-1./slope) -1 - prefactor = 2./(2*slope -3) * kappa_0 * theta_c**2 / a_factor + a_factor = 0.5**(-1./slope) - 1 + prefactor = 2./(2*slope - 3) * kappa_0 * theta_c**2 / a_factor denominator_factor = (1 + a_factor * r**2/theta_c**2)**(slope - 3./2) return prefactor/r * (1 - 1/denominator_factor) @@ -107,7 +110,7 @@ def derivatives(self, x, y, kappa_0, theta_c, center_x=0, center_y=0, slope=8): x_ = x - center_x y_ = y - center_y R = np.sqrt(x_**2 + y_**2) - R = np.maximum(R,0.00000001) + R = np.maximum(R, 0.00000001) f_x = self.alpha_radial(R, kappa_0, theta_c, slope) * x_ / R f_y = self.alpha_radial(R, kappa_0, theta_c, slope) * y_ / R return f_x, f_y @@ -126,8 +129,8 @@ def hessian(self, x, y, kappa_0, theta_c, center_x=0, center_y=0, slope=8): x_ = x - center_x y_ = y - center_y R = np.sqrt(x_**2 + y_**2) - R = np.maximum(R,0.00000001) - a_factor = (0.5)**(-1./slope) -1 + R = np.maximum(R, 0.00000001) + a_factor = 0.5**(-1./slope) - 1 prefactor = 2./(2*slope -3) * kappa_0 * theta_c**2 / a_factor # denominator factor denominator = 1 + a_factor * R**2/theta_c**2 @@ -142,6 +145,7 @@ def density(self, R, kappa_0, theta_c, slope=8): """ three dimensional ULDM profile in angular units (rho0_physical = rho0_angular Sigma_crit / D_lens) + :param R: radius of interest :param kappa_0: central convergence of profile :param theta_c: core radius (in arcsec) @@ -149,8 +153,8 @@ def density(self, R, kappa_0, theta_c, slope=8): :return: rho(R) density in angular units """ rhotilde = self.rhotilde(kappa_0, theta_c, slope) - a_factor = (0.5)**(-1./slope) -1 - return rhotilde/(1 + a_factor* (R/theta_c)**2)**slope + a_factor = 0.5**(-1./slope) - 1 + return rhotilde/(1 + a_factor * (R/theta_c)**2)**slope def density_lens(self, r, kappa_0, theta_c, slope=8): """ @@ -166,7 +170,8 @@ def density_lens(self, r, kappa_0, theta_c, slope=8): """ return self.density(r, kappa_0, theta_c, slope) - def kappa_r(self, R, kappa_0, theta_c, slope=8): + @staticmethod + def kappa_r(R, kappa_0, theta_c, slope=8): """ convergence of the cored density profile. This routine is also for testing @@ -176,16 +181,16 @@ def kappa_r(self, R, kappa_0, theta_c, slope=8): :param slope: exponent entering the profile :return: convergence at r """ - a_factor = (0.5)**(-1./slope) -1 - return kappa_0 * (1 + a_factor * (R/theta_c)**2)**(1./2 - slope) - + a_factor = (0.5)**(-1./slope) -1 + return kappa_0 * (1 + a_factor * (R/theta_c)**2)**(1./2 - slope) def density_2d(self, x, y, kappa_0, theta_c, center_x=0, center_y=0, slope=8): """ projected two dimensional ULDM profile (convergence * Sigma_crit), but given our units convention for rho0, it is basically the convergence - :param R: radius of interest + :param x: x-coordinate + :param y: y-coordinate :param kappa_0: central convergence of profile :param theta_c: core radius (in arcsec) :param slope: exponent entering the profile @@ -198,7 +203,8 @@ def density_2d(self, x, y, kappa_0, theta_c, center_x=0, center_y=0, slope=8): def _mass_integral(self, x, slope=8): """ - Returns the analitic result of the integral appearing in mass expression + Returns the analytic result of the integral appearing in mass expression + :param slope: exponent entering the profile :return: integral result """ @@ -208,6 +214,7 @@ def _mass_integral(self, x, slope=8): def mass_3d(self, R, kappa_0, theta_c, slope=8): """ mass enclosed a 3d sphere or radius r + :param R: radius in arcseconds :param kappa_0: central convergence of profile :param theta_c: core radius (in arcsec) @@ -215,7 +222,7 @@ def mass_3d(self, R, kappa_0, theta_c, slope=8): :return: mass of soliton in angular units """ rhotilde = self.rhotilde(kappa_0, theta_c, slope) - a_factor = (0.5)**(-1./slope) -1 + a_factor = 0.5**(-1./slope) - 1 prefactor = 4. * np.pi * rhotilde * theta_c**3 / (a_factor)**(1.5) m_3d = prefactor * (self._mass_integral(R/theta_c * np.sqrt(a_factor), slope) - self._mass_integral(0, slope) ) @@ -224,6 +231,7 @@ def mass_3d(self, R, kappa_0, theta_c, slope=8): def mass_3d_lens(self, r, kappa_0, theta_c, slope=8): """ mass enclosed a 3d sphere or radius r + :param r: radius over which the mass is computed :param kappa_0: central convergence of profile :param theta_c: core radius (in arcsec) @@ -236,6 +244,7 @@ def mass_3d_lens(self, r, kappa_0, theta_c, slope=8): def mass_2d(self, R, kappa_0, theta_c, slope=8): """ mass enclosed a 2d sphere or radius r + :param R: radius over which the mass is computed :param kappa_0: central convergence of profile :param theta_c: core radius (in arcsec) diff --git a/lenstronomy/LensModel/QuadOptimizer/multi_plane_fast.py b/lenstronomy/LensModel/QuadOptimizer/multi_plane_fast.py index 53b6df42f..738f41bab 100644 --- a/lenstronomy/LensModel/QuadOptimizer/multi_plane_fast.py +++ b/lenstronomy/LensModel/QuadOptimizer/multi_plane_fast.py @@ -29,7 +29,7 @@ def __init__(self, x_image, y_image, z_lens, z_source, lens_model_list, redshift :param astropy_instance: instance of astropy to pass to lens model :param param_class: an instance of ParamClass (see documentation in QuadOptimmizer.param_manager) :param foreground_rays: (optional) pre-computed foreground rays from a previous iteration, if they are not specified - they will be re-computed + they will be re-computed :param tol_source: source plane chi^2 sigma :param numerical_alpha_class: class for computing numerically tabulated deflection angles """ diff --git a/lenstronomy/LensModel/QuadOptimizer/optimizer.py b/lenstronomy/LensModel/QuadOptimizer/optimizer.py index 2b20fdc49..088f6c624 100644 --- a/lenstronomy/LensModel/QuadOptimizer/optimizer.py +++ b/lenstronomy/LensModel/QuadOptimizer/optimizer.py @@ -39,7 +39,7 @@ def __init__(self, x_image, y_image, lens_model_list, redshift_list, z_lens, z_s :param re_optimize_scale: float, controls how tight the initial spread of particles is :param pso_convergence_mean: when to terminate the PSO fit :param foreground_rays: (optional) can pass in pre-computed foreground light rays from a previous fit - so as to not waste time recomputing them + so as to not waste time recomputing them :param tol_source: sigma in the source plane chi^2 :param tol_simplex_func: tolerance for the downhill simplex optimization :param simplex_n_iterations: number of iterations per dimension for the downhill simplex optimization diff --git a/lenstronomy/LensModel/Solver/epl_shear_solver.py b/lenstronomy/LensModel/Solver/epl_shear_solver.py index 68b6762b0..73089677f 100644 --- a/lenstronomy/LensModel/Solver/epl_shear_solver.py +++ b/lenstronomy/LensModel/Solver/epl_shear_solver.py @@ -1,4 +1,4 @@ -__author__ = 'ewoudwempe' +__author__ = 'ewoudwempe', 'sibirrer' import numpy as np from lenstronomy.LensModel.Util.epl_util import min_approx, pol_to_cart, cart_to_pol, cdot, ps, rotmat, solvequadeq, brentq_inline @@ -7,6 +7,7 @@ from lenstronomy.LensModel.Profiles.epl_numba import alpha, omega from lenstronomy.Util.numba_util import jit from lenstronomy.Util.param_util import ellipticity2phi_q, shear_cartesian2polar, shear_polar2cartesian +from lenstronomy.LensModel.Profiles.shear import Shear @jit() @@ -88,6 +89,7 @@ def _getphi(thpl, args): """ Finds all roots to both versions the 1-dimensional lens equation in phi, by doing a grid search for sign changes on the supplied thpl. In the case of extrema, refine at the relevant location. + :param thpl: What points to calculate the equation on use for detecting sign changes :param args: Parameters to be passed to the lens equation :return: an array containing all roots @@ -160,9 +162,16 @@ def solvelenseq_majoraxis(args, Nmeas=200, Nmeas_extra=50): def _check_center(kwargs_lens): """Checks if the shear-at-center convention is properly used.""" - if kwargs_lens[1]['ra_0'] != kwargs_lens[0]['center_x'] or kwargs_lens[1]['dec_0'] != kwargs_lens[0]['center_y']: - raise ValueError("Center of lens (center_{x,y}) must be the same as center of shear ({ra,dec}_0). " - "This can be ensured by supplying a dictionary-style joint_setting_list to the model.") + # calculate (inverse) displacement caused by the offset between shear and lens centroid + # this shift needs to be added to the source position such that the solution of the lens equation + # without this shift in the shear is the correct one + if len(kwargs_lens) > 1: + shear = Shear() + # calculate shift from the deflector centroid from the shear field + alpha_x, alpha_y = shear.derivatives(kwargs_lens[0]['center_x'], kwargs_lens[0]['center_y'], **kwargs_lens[1]) + return alpha_x, alpha_y + else: + return 0, 0 def solve_lenseq_pemd(pos_, kwargs_lens, Nmeas=400, Nmeas_extra=80, **kwargs): @@ -182,14 +191,15 @@ def solve_lenseq_pemd(pos_, kwargs_lens, Nmeas=400, Nmeas_extra=80, **kwargs): theta_ell, q = ellipticity2phi_q(kwargs_lens[0]['e1'], kwargs_lens[0]['e2']) b = kwargs_lens[0]['theta_E']*np.sqrt(q) - - cen = kwargs_lens[0]['center_x']+1j*kwargs_lens[0]['center_y'] - p = pos[0]+1j*pos[1]-cen if len(kwargs_lens) > 1: gamma = kwargs_lens[1]['gamma1']+1j*kwargs_lens[1]['gamma2'] - _check_center(kwargs_lens) else: gamma = 0+0j + shift_x, shift_y = _check_center(kwargs_lens) + shift = shift_x + 1j * shift_y + cen = kwargs_lens[0]['center_x']+1j*kwargs_lens[0]['center_y'] + p = pos[0]+1j*pos[1] - cen + shift + rotfact = np.exp(-1j*theta_ell) gamma *= rotfact**2 p *= rotfact @@ -218,7 +228,6 @@ def caustics_epl_shear(kwargs_lens, num_th=500, maginf=0, sourceplane=True, retu gamma1unr, gamma2unr = kwargs_lens[1]['gamma1'], kwargs_lens[1]['gamma2'] else: gamma1unr, gamma2unr = 0, 0 - _check_center(kwargs_lens) t = kwargs_lens[0]['gamma']-1 if 'gamma' in kwargs_lens[0] else 1 theta_ell, q = ellipticity2phi_q(e1, e2) theta_gamma, gamma_mag = shear_cartesian2polar(gamma1unr, gamma2unr) diff --git a/lenstronomy/LensModel/Solver/lens_equation_solver.py b/lenstronomy/LensModel/Solver/lens_equation_solver.py index 3782b3588..64f0931cd 100644 --- a/lenstronomy/LensModel/Solver/lens_equation_solver.py +++ b/lenstronomy/LensModel/Solver/lens_equation_solver.py @@ -100,7 +100,7 @@ def candidate_solutions(self, sourcePos_x, sourcePos_y, kwargs_lens, min_distanc absmapped = util.displaceAbs(x_mapped, y_mapped, sourcePos_x, sourcePos_y) # select minima in the grid points and select grid points that do not deviate more than the # width of the grid point to a solution of the lens equation - x_mins, y_mins, delta_map = util.neighborSelect(absmapped, x_grid, y_grid) + x_mins, y_mins, delta_map = util.local_minima_2d(absmapped, x_grid, y_grid) # pixel width pixel_width = x_grid[1]-x_grid[0] @@ -124,7 +124,7 @@ def image_position_analytical(self, x, y, kwargs_lens, arrival_time_sort=True, m """ lens_model_list = list(self.lensModel.lens_model_list) if lens_model_list not in (['SIE', 'SHEAR'], ['SIE'], ['EPL_NUMBA', 'SHEAR'], ['EPL_NUMBA'], ['EPL', 'SHEAR'], ['EPL']): - raise ValueError("Only SIE or PEMD (+shear) supported in the analytical solver for now") + raise ValueError("Only SIE, EPL, EPL_NUMBA (+shear) supported in the analytical solver for now.") x_mins, y_mins = solve_lenseq_pemd((x, y), kwargs_lens, **kwargs_solver) if arrival_time_sort: @@ -189,7 +189,8 @@ def image_position_lenstronomy(self, sourcePos_x, sourcePos_y, kwargs_lens, min_ # pixel width x_mins, y_mins, delta_map, pixel_width = self.candidate_solutions(sourcePos_x, sourcePos_y, kwargs_lens, min_distance, search_window, verbose, x_center, y_center) if verbose: - print("There are %s regions identified that could contain a solution of the lens equation" % len(x_mins)) + print("There are %s regions identified that could contain a solution of the lens equation with" + "coordinates %s and %s " % (len(x_mins), x_mins, y_mins)) if len(x_mins) < 1: return x_mins, y_mins if initial_guess_cut: diff --git a/lenstronomy/LensModel/lens_model.py b/lenstronomy/LensModel/lens_model.py index 2c603057c..f4c82d773 100644 --- a/lenstronomy/LensModel/lens_model.py +++ b/lenstronomy/LensModel/lens_model.py @@ -22,28 +22,28 @@ def __init__(self, lens_model_list, z_lens=None, z_source=None, lens_redshift_li :param lens_model_list: list of strings with lens model names :param z_lens: redshift of the deflector (only considered when operating in single plane mode). - Is only needed for specific functions that require a cosmology. + Is only needed for specific functions that require a cosmology. :param z_source: redshift of the source: Needed in multi_plane option only, - not required for the core functionalities in the single plane mode. + not required for the core functionalities in the single plane mode. :param lens_redshift_list: list of deflector redshift (corresponding to the lens model list), - only applicable in multi_plane mode. + only applicable in multi_plane mode. :param cosmo: instance of the astropy cosmology class. If not specified, uses the default cosmology. :param multi_plane: bool, if True, uses multi-plane mode. Default is False. :param numerical_alpha_class: an instance of a custom class for use in NumericalAlpha() lens model - (see documentation in Profiles/numerical_alpha) + (see documentation in Profiles/numerical_alpha) :param kwargs_interp: interpolation keyword arguments specifying the numerics. See description in the Interpolate() class. Only applicable for 'INTERPOL' and 'INTERPOL_SCALED' models. :param observed_convention_index: a list of indices, corresponding to the lens_model_list element with same - index, where the 'center_x' and 'center_y' kwargs correspond to observed (lensed) positions, not physical - positions. The code will compute the physical locations when performing computations + index, where the 'center_x' and 'center_y' kwargs correspond to observed (lensed) positions, not physical + positions. The code will compute the physical locations when performing computations :param z_source_convention: float, redshift of a source to define the reduced deflection angles of the lens - models. If None, 'z_source' is used. + models. If None, 'z_source' is used. :param cosmo_interp: boolean (only employed in multi-plane mode), interpolates astropy.cosmology distances for - faster calls when accessing several lensing planes + faster calls when accessing several lensing planes :param z_interp_stop: (only in multi-plane with cosmo_interp=True); maximum redshift for distance interpolation - This number should be higher or equal the maximum of the source redshift and/or the z_source_convention + This number should be higher or equal the maximum of the source redshift and/or the z_source_convention :param num_z_interp: (only in multi-plane with cosmo_interp=True); number of redshift bins for interpolating - distances + distances """ self.lens_model_list = lens_model_list self.z_lens = z_lens diff --git a/lenstronomy/LensModel/lens_model_extensions.py b/lenstronomy/LensModel/lens_model_extensions.py index f829d7ed4..d201d2192 100644 --- a/lenstronomy/LensModel/lens_model_extensions.py +++ b/lenstronomy/LensModel/lens_model_extensions.py @@ -13,12 +13,12 @@ def __init__(self, lensModel): """ :param lensModel: instance of the LensModel() class, or with same functionalities. - In particular, the following definitions are required to execute all functionalities presented in this class: - def ray_shooting() - def magnification() - def kappa() - def alpha() - def hessian() + In particular, the following definitions are required to execute all functionalities presented in this class: + def ray_shooting() + def magnification() + def kappa() + def alpha() + def hessian() """ self._lensModel = lensModel diff --git a/lenstronomy/LensModel/profile_list_base.py b/lenstronomy/LensModel/profile_list_base.py index c1570729b..569427546 100644 --- a/lenstronomy/LensModel/profile_list_base.py +++ b/lenstronomy/LensModel/profile_list_base.py @@ -31,7 +31,7 @@ def __init__(self, lens_model_list, numerical_alpha_class=None, lens_redshift_li :param lens_model_list: list of strings with lens model names :param numerical_alpha_class: an instance of a custom class for use in NumericalAlpha() lens model - deflection angles as a lens model. See the documentation in Profiles.numerical_deflections + deflection angles as a lens model. See the documentation in Profiles.numerical_deflections :param kwargs_interp: interpolation keyword arguments specifying the numerics. See description in the Interpolate() class. Only applicable for 'INTERPOL' and 'INTERPOL_SCALED' models. """ diff --git a/lenstronomy/LightModel/Profiles/interpolation.py b/lenstronomy/LightModel/Profiles/interpolation.py index aed35addb..2c438aba1 100644 --- a/lenstronomy/LightModel/Profiles/interpolation.py +++ b/lenstronomy/LightModel/Profiles/interpolation.py @@ -13,7 +13,7 @@ class Interpol(object): class which uses an interpolation of an image to compute the surface brightness parameters are - 'image': 2d numpy array of surface brightness (not integrated flux per pixel!) + 'image': 2d numpy array of surface brightness per square arc second (not integrated flux per pixel!) 'center_x': coordinate of center of image in angular units (i.e. arc seconds) 'center_y': coordinate of center of image in angular units (i.e. arc seconds) 'phi_G': rotation of image relative to the rectangular ra-to-dec orientation @@ -32,8 +32,10 @@ def function(self, x, y, image=None, amp=1, center_x=0, center_y=0, phi_G=0, sca :param x: x-coordinate to evaluate surface brightness :param y: y-coordinate to evaluate surface brightness - :param image: 2d numpy array (image) to be used to interpolate - :param amp: amplitude of surface brightness scaling in respect of original image + :param image: pixelized surface brightness (an image) to be used to interpolate in units of surface brightness + (flux per square arc seconds, not flux per pixel!) + :type image: 2d numpy array + :param amp: amplitude of surface brightness scaling in respect of original input image :param center_x: center of interpolated image :param center_y: center of interpolated image :param phi_G: rotation angle of simulated image in respect to input gird @@ -60,12 +62,14 @@ def image_interp(self, x, y, image): # (try reversing, the unit tests will fail) return self._image_interp(y, x, grid=False) - def total_flux(self, image, scale, amp=1, center_x=0, center_y=0, phi_G=0): + @staticmethod + def total_flux(image, scale, amp=1, center_x=0, center_y=0, phi_G=0): """ - sums up all the image surface brightness (image pixels defined in surface brightness at the coordinate of the pixel) - times pixel area + sums up all the image surface brightness (image pixels defined in surface brightness at the coordinate of the + pixel) times pixel area - :param image: pixelized surface brightness + :param image: pixelized surface brightness used to interpolate in units of surface brightness + (flux per square arc seconds, not flux per pixel!) :param scale: scale of the pixel in units of angle :param amp: linear scaling parameter of the surface brightness multiplicative with the initial image :param center_x: center of image in angular coordinates diff --git a/lenstronomy/LightModel/Profiles/sersic.py b/lenstronomy/LightModel/Profiles/sersic.py index 5f0cabeea..cb77fdfd7 100644 --- a/lenstronomy/LightModel/Profiles/sersic.py +++ b/lenstronomy/LightModel/Profiles/sersic.py @@ -51,6 +51,7 @@ class SersicElliptic(SersicUtil): this class contains functions to evaluate an elliptical Sersic function .. math:: + I(R) = I_0 \\exp \\left[ -b_n (R/R_{\\rm Sersic})^{\\frac{1}{n}}\\right] with :math:`I_0 = amp`, @@ -93,8 +94,9 @@ class CoreSersic(SersicUtil): this class contains the Core-Sersic function introduced by e.g Trujillo et al. 2004 .. math:: + I(R) = I' \\left[1 + (R_b/R)^{\\alpha} \\right]^{\\gamma / \\alpha} - \\exp \\left{ -b_n \\left[(R^{\\alpha} + R_b^{\alpha})/R_e^{\\alpha} \\right]^{1 / (n\\alpha)} \\right} + \\exp \\left{ -b_n \\left[(R^{\\alpha} + R_b^{\\alpha})/R_e^{\\alpha} \\right]^{1 / (n\\alpha)} \\right} with diff --git a/lenstronomy/LightModel/Profiles/shapelets.py b/lenstronomy/LightModel/Profiles/shapelets.py index b96688cd1..96d76341e 100644 --- a/lenstronomy/LightModel/Profiles/shapelets.py +++ b/lenstronomy/LightModel/Profiles/shapelets.py @@ -284,6 +284,7 @@ def shapelet_basis_2d(self, num_order, beta, numPix, deltaPix=1, center_x=0, cen def decomposition(self, image, x, y, n_max, beta, deltaPix, center_x=0, center_y=0): """ decomposes an image into the shapelet coefficients in same order as for the function call + :param image: :param x: :param y: diff --git a/lenstronomy/LightModel/Profiles/starlets_util.py b/lenstronomy/LightModel/Profiles/starlets_util.py index 0731886e8..89c647cce 100644 --- a/lenstronomy/LightModel/Profiles/starlets_util.py +++ b/lenstronomy/LightModel/Profiles/starlets_util.py @@ -1,7 +1,7 @@ -__author__ = 'herjy', 'aymgal' +__author__ = 'herjy', 'aymgal', 'sibirrer' import numpy as np -import scipy.ndimage.filters as scf +from scipy import ndimage from lenstronomy.Util.package_util import exporter export, __all__ = exporter() @@ -47,17 +47,17 @@ def transform(img, n_scales, second_gen=False): ######Calculates c(j+1) ###### Line convolution - cnew = scf.convolve1d(c, newh[0, :], axis=0, mode=mode) + cnew = ndimage.convolve1d(c, newh[0, :], axis=0, mode=mode) ###### Column convolution - cnew = scf.convolve1d(cnew, newh[0, :], axis=1, mode=mode) + cnew = ndimage.convolve1d(cnew, newh[0, :], axis=1, mode=mode) if second_gen: ###### hoh for g; Column convolution - hc = scf.convolve1d(cnew, newh[0, :], axis=0, mode=mode) + hc = ndimage.convolve1d(cnew, newh[0, :], axis=0, mode=mode) ###### hoh for g; Line convolution - hc = scf.convolve1d(hc, newh[0, :], axis=1, mode=mode) + hc = ndimage.convolve1d(hc, newh[0, :], axis=1, mode=mode) ###### wj+1 = cj - hcj+1 wave[i, :, :] = c - hc @@ -101,9 +101,9 @@ def inverse_transform(wave, fast=True, second_gen=False): H = np.dot(newh.T, newh) ###### Line convolution - cnew = scf.convolve1d(cJ, newh[0, :], axis=0, mode=mode) + cnew = ndimage.convolve1d(cJ, newh[0, :], axis=0, mode=mode) ###### Column convolution - cnew = scf.convolve1d(cnew, newh[0, :], axis=1, mode=mode) + cnew = ndimage.convolve1d(cnew, newh[0, :], axis=1, mode=mode) cJ = cnew + wave[lvl-1-i, :, :] diff --git a/lenstronomy/LightModel/light_model.py b/lenstronomy/LightModel/light_model.py index 2babf2538..07ad0cd2c 100644 --- a/lenstronomy/LightModel/light_model.py +++ b/lenstronomy/LightModel/light_model.py @@ -28,7 +28,7 @@ def __init__(self, light_model_list, deflection_scaling_list=None, source_redshi :param light_model_list: list of light models :param deflection_scaling_list: list of floats indicating a relative scaling of the deflection angle from the reduced angles in the lens model definition (optional, only possible in single lens plane with multiple source - planes) + planes) :param source_redshift_list: list of redshifts for the different light models (optional and only used in multi-plane lensing in conjunction with a cosmology model) :param smoothing: smoothing factor for certain models (deprecated) diff --git a/lenstronomy/Plots/plot_quasar_images.py b/lenstronomy/Plots/plot_quasar_images.py index a4ae9c236..a8f034e80 100644 --- a/lenstronomy/Plots/plot_quasar_images.py +++ b/lenstronomy/Plots/plot_quasar_images.py @@ -23,19 +23,19 @@ def plot_quasar_images(lens_model, x_image, y_image, source_x, source_y, kwargs_ :param z_source: the source redshift :param cosmo: (optional) an instance of astropy.cosmology; if not specified, a default cosmology will be used :param grid_resolution: the grid resolution in units arcsec/pixel; if not specified, an appropriate value will - be estimated from the source size + be estimated from the source size :param grid_radius_arcsec: (optional) the size of the ray tracing region in arcsec; if not specified, an appropriate value - will be estimated from the source size + will be estimated from the source size :param source_light_model: the model for background source light; currently implemented are 'SINGLE_GAUSSIAN' and - 'DOUBLE_GAUSSIAN'. + 'DOUBLE_GAUSSIAN'. :param dx: used with source model 'DOUBLE_GAUSSIAN', the offset of the second source light profile from the first - [arcsec] + [arcsec] :param dy: used with source model 'DOUBLE_GAUSSIAN', the offset of the second source light profile from the first - [arcsec] + [arcsec] :param size_scale: used with source model 'DOUBLE_GAUSSIAN', the size of the second source light profile relative - to the first + to the first :param amp_scale: used with source model 'DOUBLE_GAUSSIAN', the peak brightness of the second source light profile - relative to the first + relative to the first :return: Four images of the background source in the image plane """ diff --git a/lenstronomy/Plots/plot_util.py b/lenstronomy/Plots/plot_util.py index 9e796ce2d..f830ff890 100644 --- a/lenstronomy/Plots/plot_util.py +++ b/lenstronomy/Plots/plot_util.py @@ -10,14 +10,14 @@ def sqrt(inputArray, scale_min=None, scale_max=None): """Performs sqrt scaling of the input numpy array. - @type inputArray: numpy array - @param inputArray: image data array - @type scale_min: float - @param scale_min: minimum data value - @type scale_max: float - @param scale_max: maximum data value - @rtype: numpy array - @return: image data array + :type inputArray: numpy array + :param inputArray: image data array + :type scale_min: float + :param scale_min: minimum data value + :type scale_max: float + :param scale_max: maximum data value + :rtype: numpy array + :return: image data array """ diff --git a/lenstronomy/PointSource/point_source.py b/lenstronomy/PointSource/point_source.py index 7e3169495..e59308574 100644 --- a/lenstronomy/PointSource/point_source.py +++ b/lenstronomy/PointSource/point_source.py @@ -21,7 +21,7 @@ def __init__(self, point_source_type_list, lensModel=None, fixed_magnification_l This option then requires to provide a 'source_amp' amplitude of the source brightness instead of 'point_amp' the list of image brightnesses. :param additional_images_list: list of booleans (same length as point_source_type_list). If True, search for - additional images of the same source is conducted. + additional images of the same source is conducted. :param flux_from_point_source_list: list of booleans (optional), if set, will only return image positions (for imaging modeling) for the subset of the point source lists that =True. This option enables to model imaging data with transient point sources, when the point source positions are measured and present at a @@ -32,10 +32,8 @@ def __init__(self, point_source_type_list, lensModel=None, fixed_magnification_l equation is conducted with the lens model parameters provided. This can increase the speed as multiple times the image positions are requested for the same lens model. Attention in usage! :param kwargs_lens_eqn_solver: keyword arguments specifying the numerical settings for the lens equation solver - see LensEquationSolver() class for details - - for the kwargs_lens_eqn_solver parameters: have a look at the lensEquationSolver class, such as: - min_distance=0.01, search_window=5, precision_limit=10**(-10), num_iter_max=100 + see LensEquationSolver() class for details ,such as: + min_distance=0.01, search_window=5, precision_limit=10**(-10), num_iter_max=100 """ self._lensModel = lensModel diff --git a/lenstronomy/PointSource/point_source_param.py b/lenstronomy/PointSource/point_source_param.py index 890deb944..84bf989e5 100644 --- a/lenstronomy/PointSource/point_source_param.py +++ b/lenstronomy/PointSource/point_source_param.py @@ -2,10 +2,59 @@ __all__ = ['PointSourceParam'] +from lenstronomy.Sampling.param_group import ModelParamGroup, SingleParam, ArrayParam -class PointSourceParam(object): + +class SourcePositionParam(SingleParam): + """ + Source position parameter, ra_source and dec_source """ + param_names = ['ra_source', 'dec_source'] + _kwargs_lower = {'ra_source': -100, 'dec_source': -100} + _kwargs_upper = {'ra_source': 100, 'dec_source': 100} + +class LensedPosition(ArrayParam): + """ + Represents lensed positions, possibly many. ra_image and dec_image + + :param num_images: integer. The number of lensed positions to model. + """ + _kwargs_lower = {'ra_image': -100, 'dec_image': -100, } + _kwargs_upper = {'ra_image': 100, 'dec_image': 100, } + + def __init__(self, num_images): + ArrayParam.__init__(self, int(num_images) > 0) + self.param_names = {'ra_image': int(num_images), 'dec_image': int(num_images)} + + +class SourceAmp(SingleParam): + """ + Source amplification + """ + param_names = ['source_amp'] + _kwargs_lower = {'source_amp': 0} + _kwargs_upper = {'source_amp': 100} + + +class ImageAmp(ArrayParam): + """ + Observed amplification of lensed images of a point source. Can model + arbitrarily many magnified images + + :param num_point_sources: integer. The number of lensed images without fixed magnification. + """ + _kwargs_lower = {'point_amp': 0} + _kwargs_upper = {'point_amp': 100} + + def __init__(self, num_point_sources): + ArrayParam.__init__(self, int(num_point_sources) > 0) + self.param_names = {'point_amp': int(num_point_sources)} + + +class PointSourceParam(object): + """ + Point source parameters """ def __init__(self, model_list, kwargs_fixed, num_point_source_list=None, linear_solver=True, @@ -16,9 +65,9 @@ def __init__(self, model_list, kwargs_fixed, num_point_source_list=None, linear_ :param kwargs_fixed: list of keyword arguments with parameters to be held fixed :param num_point_source_list: list of number of point sources per point source model class :param linear_solver: bool, if True, does not return linear parameters for the sampler - (will be solved linearly instead) + (will be solved linearly instead) :param fixed_magnification_list: list of booleans, if entry is True, keeps one overall scaling among the - point sources in this class + point sources in this class """ self.model_list = model_list if num_point_source_list is None: @@ -32,36 +81,40 @@ def __init__(self, model_list, kwargs_fixed, num_point_source_list=None, linear_ self.kwargs_fixed = self.add_fix_linear(kwargs_fixed) self._linear_solver = linear_solver + self.param_groups = [] + for i, model in enumerate(self.model_list): + params = [] + num = num_point_source_list[i] + if model in ['LENSED_POSITION', 'UNLENSED']: + params.append(LensedPosition(num)) + elif model == 'SOURCE_POSITION': + params.append(SourcePositionParam(True)) + else: + raise ValueError("%s not a valid point source model" % model) + + if fixed_magnification_list[i] and model in ['LENSED_POSITION', 'SOURCE_POSITION']: + params.append(SourceAmp(True)) + else: + params.append(ImageAmp(num)) + + self.param_groups.append(params) + if kwargs_lower is None: kwargs_lower = [] - for k, model in enumerate(self.model_list): - num = self._num_point_sources_list[k] - if model in ['LENSED_POSITION', 'UNLENSED']: - fixed_low = {'ra_image': [-100] * num, 'dec_image': [-100] * num} - elif model in ['SOURCE_POSITION']: - fixed_low = {'ra_source': -100, 'dec_source': -100} - else: - raise ValueError("%s not a valid point source model" % model) - if self._fixed_magnification_list[k] is True and model in ['LENSED_POSITION', 'SOURCE_POSITION']: - fixed_low['source_amp'] = 0 - else: - fixed_low['point_amp'] = np.zeros(num) - kwargs_lower.append(fixed_low) + for model_params in self.param_groups: + fixed_lower = {} + for param_group in model_params: + fixed_lower = dict(fixed_lower, **param_group.kwargs_lower) + kwargs_lower.append(fixed_lower) + if kwargs_upper is None: kwargs_upper = [] - for k, model in enumerate(self.model_list): - num = self._num_point_sources_list[k] - if model in ['LENSED_POSITION', 'UNLENSED']: - fixed_high = {'ra_image': [100] * num, 'dec_image': [100] * num} - elif model in ['SOURCE_POSITION']: - fixed_high = {'ra_source': 100, 'dec_source': 100} - else: - raise ValueError("%s not a valid point source model" % model) - if self._fixed_magnification_list[k] is True and model in ['LENSED_POSITION', 'SOURCE_POSITION']: - fixed_high['source_amp'] = 100 - else: - fixed_high['point_amp'] = np.ones(num)*100 - kwargs_upper.append(fixed_high) + for model_params in self.param_groups: + fixed_upper = {} + for param_group in model_params: + fixed_upper = dict(fixed_upper, **param_group.kwargs_upper) + kwargs_upper.append(fixed_upper) + self.lower_limit = kwargs_lower self.upper_limit = kwargs_upper @@ -73,45 +126,10 @@ def get_params(self, args, i): :return: keyword argument list of point sources, index relevant for the next class """ kwargs_list = [] - for k, model in enumerate(self.model_list): - kwargs = {} - kwargs_fixed = self.kwargs_fixed[k] - if model in ['LENSED_POSITION', 'UNLENSED']: - if 'ra_image' not in kwargs_fixed: - kwargs['ra_image'] = np.array(args[i:i + self._num_point_sources_list[k]]) - i += self._num_point_sources_list[k] - else: - kwargs['ra_image'] = kwargs_fixed['ra_image'] - if 'dec_image' not in kwargs_fixed: - kwargs['dec_image'] = np.array(args[i:i + self._num_point_sources_list[k]]) - i += self._num_point_sources_list[k] - else: - kwargs['dec_image'] = kwargs_fixed['dec_image'] - if model in ['SOURCE_POSITION']: - if 'ra_source' not in kwargs_fixed: - kwargs['ra_source'] = args[i] - i += 1 - else: - kwargs['ra_source'] = kwargs_fixed['ra_source'] - if 'dec_source' not in kwargs_fixed: - kwargs['dec_source'] = args[i] - i += 1 - else: - kwargs['dec_source'] = kwargs_fixed['dec_source'] - # amplitude parameter handling - if self._fixed_magnification_list[k] is True and model in ['LENSED_POSITION', 'SOURCE_POSITION']: - if 'source_amp' not in kwargs_fixed: - kwargs['source_amp'] = args[i] - i += 1 - else: - kwargs['source_amp'] = kwargs_fixed['source_amp'] - else: - if 'point_amp' not in kwargs_fixed: - kwargs['point_amp'] = np.array(args[i:i + self._num_point_sources_list[k]]) - i += self._num_point_sources_list[k] - else: - kwargs['point_amp'] = kwargs_fixed['point_amp'] - + for k, param_group in enumerate(self.param_groups): + kwargs, i = ModelParamGroup.compose_get_params( + param_group, args, i, kwargs_fixed=self.kwargs_fixed[k] + ) kwargs_list.append(kwargs) return kwargs_list, i @@ -122,32 +140,12 @@ def set_params(self, kwargs_list): :return: sorted list of parameters being sampled extracted from kwargs_list """ args = [] - for k, model in enumerate(self.model_list): + for k, param_group in enumerate(self.param_groups): kwargs = kwargs_list[k] kwargs_fixed = self.kwargs_fixed[k] - if model in ['LENSED_POSITION', 'UNLENSED']: - if 'ra_image' not in kwargs_fixed: - x_pos = kwargs['ra_image'][0:self._num_point_sources_list[k]] - for x in x_pos: - args.append(x) - if 'dec_image' not in kwargs_fixed: - y_pos = kwargs['dec_image'][0:self._num_point_sources_list[k]] - for y in y_pos: - args.append(y) - if model in ['SOURCE_POSITION']: - if 'ra_source' not in kwargs_fixed: - args.append(kwargs['ra_source']) - if 'dec_source' not in kwargs_fixed: - args.append(kwargs['dec_source']) - # amplitude parameter handling - if self._fixed_magnification_list[k] is True and model in ['LENSED_POSITION', 'SOURCE_POSITION']: - if 'source_amp' not in kwargs_fixed: - args.append(kwargs['source_amp']) - else: - if 'point_amp' not in kwargs_fixed: - amp = kwargs['point_amp'][0:self._num_point_sources_list[k]] - for a in amp: - args.append(a) + args.extend(ModelParamGroup.compose_set_params( + param_group, kwargs, kwargs_fixed=kwargs_fixed + )) return args def num_param(self): @@ -156,36 +154,13 @@ def num_param(self): :return: int, list of parameter names """ - num = 0 - name_list = [] - for k, model in enumerate(self.model_list): - kwargs_fixed = self.kwargs_fixed[k] - if model in ['LENSED_POSITION', 'UNLENSED']: - if 'ra_image' not in kwargs_fixed: - num += self._num_point_sources_list[k] - for i in range(self._num_point_sources_list[k]): - name_list.append('ra_image') - if 'dec_image' not in kwargs_fixed: - num += self._num_point_sources_list[k] - for i in range(self._num_point_sources_list[k]): - name_list.append('dec_image') - if model in ['SOURCE_POSITION']: - if 'ra_source' not in kwargs_fixed: - num += 1 - name_list.append('ra_source') - if 'dec_source' not in kwargs_fixed: - num += 1 - name_list.append('dec_source') - # amplitude handling - if self._fixed_magnification_list[k] is True and model in ['LENSED_POSITION', 'SOURCE_POSITION']: - if 'source_amp' not in kwargs_fixed: - num += 1 - name_list.append('source_amp') - else: - if 'point_amp' not in kwargs_fixed: - num += self._num_point_sources_list[k] - for i in range(self._num_point_sources_list[k]): - name_list.append('point_amp') + num, name_list = 0, [] + for k, param_group in enumerate(self.param_groups): + n, names = ModelParamGroup.compose_num_params( + param_group, kwargs_fixed=self.kwargs_fixed[k] + ) + num += n + name_list += names return num, name_list def add_fix_linear(self, kwargs_fixed): diff --git a/lenstronomy/Sampling/Likelihoods/image_likelihood.py b/lenstronomy/Sampling/Likelihoods/image_likelihood.py index 6dc48e758..fb7734a77 100644 --- a/lenstronomy/Sampling/Likelihoods/image_likelihood.py +++ b/lenstronomy/Sampling/Likelihoods/image_likelihood.py @@ -19,9 +19,9 @@ def __init__(self, multi_band_list, multi_band_type, kwargs_model, bands_compute :param image_likelihood_mask_list: list of boolean 2d arrays of size of images marking the pixels to be evaluated in the likelihood :param source_marg: marginalization addition on the imaging likelihood based on the covariance of the inferred - linear coefficients + linear coefficients :param linear_prior: float or list of floats (when multi-linear setting is chosen) indicating the range of - linear amplitude priors when computing the marginalization term. + linear amplitude priors when computing the marginalization term. :param check_positive_flux: bool, option to punish models that do not have all positive linear amplitude parameters :param kwargs_pixelbased: keyword arguments with various settings related to the pixel-based solver diff --git a/lenstronomy/Sampling/Likelihoods/position_likelihood.py b/lenstronomy/Sampling/Likelihoods/position_likelihood.py index cedbf8798..f5baefc86 100644 --- a/lenstronomy/Sampling/Likelihoods/position_likelihood.py +++ b/lenstronomy/Sampling/Likelihoods/position_likelihood.py @@ -19,7 +19,7 @@ def __init__(self, point_source_class, image_position_uncertainty=0.005, astrome :param image_position_uncertainty: uncertainty in image position uncertainty (1-sigma Gaussian radially), this is applicable for astrometric uncertainties as well as if image positions are provided as data :param astrometric_likelihood: bool, if True, evaluates the astrometric uncertainty of the predicted and modeled - image positions with an offset 'delta_x_image' and 'delta_y_image' + image positions with an offset 'delta_x_image' and 'delta_y_image' :param image_position_likelihood: bool, if True, evaluates the likelihood of the model predicted image position given the data/measured image positions :param ra_image_list: list or RA image positions per model component @@ -32,11 +32,11 @@ def __init__(self, point_source_class, image_position_uncertainty=0.005, astrome :param source_position_sigma: r.m.s. value corresponding to a 1-sigma Gaussian likelihood accepted by the model precision in matching the source position :param force_no_add_image: bool, if True, will punish additional images appearing in the frame of the modelled - image(first calculate them) + image(first calculate them) :param restrict_image_number: bool, if True, searches for all appearing images in the frame of the data and - compares with max_num_images + compares with max_num_images :param max_num_images: integer, maximum number of appearing images. Default is the number of images given in - the Param() class + the Param() class """ self._pointSource = point_source_class # TODO replace with public function of ray_shooting diff --git a/lenstronomy/Sampling/Likelihoods/prior_likelihood.py b/lenstronomy/Sampling/Likelihoods/prior_likelihood.py index 3e456ded7..b1383f90f 100644 --- a/lenstronomy/Sampling/Likelihoods/prior_likelihood.py +++ b/lenstronomy/Sampling/Likelihoods/prior_likelihood.py @@ -34,8 +34,7 @@ def __init__(self, prior_lens=None, prior_source=None, prior_lens_light=None, pr :param prior_special_kde: list of [param_name, samples] :param prior_extinction_kde: list of [index_model, param_name, samples] - :param prior_lens_lognormal: list of [index_model, param_name, mean, 1-sigma - priors] + :param prior_lens_lognormal: list of [index_model, param_name, mean, 1-sigma priors] :param prior_source_lognormal: list of [index_model, param_name, mean, 1-sigma priors] :param prior_lens_light_lognormal: list of [index_model, param_name, mean, 1-sigma priors] :param prior_ps_lognormal: list of [index_model, param_name, mean, 1-sigma priors] diff --git a/lenstronomy/Sampling/Likelihoods/time_delay_likelihood.py b/lenstronomy/Sampling/Likelihoods/time_delay_likelihood.py index a302fdd41..598e7475a 100644 --- a/lenstronomy/Sampling/Likelihoods/time_delay_likelihood.py +++ b/lenstronomy/Sampling/Likelihoods/time_delay_likelihood.py @@ -13,10 +13,10 @@ def __init__(self, time_delays_measured, time_delays_uncertainties, lens_model_c :param time_delays_measured: relative time delays (in days) in respect to the first image of the point source :param time_delays_uncertainties: time-delay uncertainties in same order as time_delay_measured. Alternatively - a full covariance matrix that describes the likelihood. + a full covariance matrix that describes the likelihood. :param lens_model_class: instance of the LensModel() class :param point_source_class: instance of the PointSource() class, note: the first point source type is the one the - time delays are imposed on + time delays are imposed on """ if time_delays_measured is None: diff --git a/lenstronomy/Sampling/Pool/multiprocessing.py b/lenstronomy/Sampling/Pool/multiprocessing.py index 66de7dff9..e28d0f8d4 100644 --- a/lenstronomy/Sampling/Pool/multiprocessing.py +++ b/lenstronomy/Sampling/Pool/multiprocessing.py @@ -43,25 +43,20 @@ class MultiPool(Pool): A modified version of :class:`multiprocessing.pool.Pool` that has better behavior with regard to ``KeyboardInterrupts`` in the :func:`map` method. (Original author: `Peter K. G. Williams `_) - - Parameters - ---------- - processes : int, optional - The number of worker processes to use; defaults to the number of CPUs. - - initializer : callable, optional - If specified, a callable that will be invoked by each worker process when it starts. - - initargs : iterable, optional - Arguments for ``initializer``; it will be called as ``initializer(*initargs)``. - - kwargs: - Extra arguments passed to the :class:`multiprocessing.pool.Pool` superclass. - """ wait_timeout = 3600 def __init__(self, processes=None, initializer=None, initargs=(), **kwargs): + """ + + :param processes: The number of worker processes to use; defaults to the number of CPUs. + :type processes: int, optional + :param initializer: If specified, a callable that will be invoked by each worker process when it starts. + :type initializer: callable, optional + :param initargs: Arguments for ``initializer``; it will be called as ``initializer(*initargs)``. + :type initargs: iterable, optional + :param kwargs: Extra arguments passed to the :class:`multiprocessing.pool.Pool` superclass. + """ new_initializer = functools.partial(_initializer_wrapper, initializer) super(MultiPool, self).__init__(processes, new_initializer, initargs, **kwargs) @@ -84,28 +79,22 @@ def map(self, func, iterable, chunksize=None, callback=None): :meth:`multiprocessing.pool.Pool.map()`, without catching ``KeyboardInterrupt``. - Parameters - ---------- - func : callable - A function or callable object that is executed on each element of + :param func: A function or callable object that is executed on each element of the specified ``tasks`` iterable. This object must be picklable (i.e. it can't be a function scoped within a function or a ``lambda`` function). This should accept a single positional argument and return a single object. - iterable : iterable - A list or iterable of tasks. Each task can be itself an iterable + :type func: callable + :param iterable: A list or iterable of tasks. Each task can be itself an iterable (e.g., tuple) of values or data to pass in to the worker function. - callback : callable, optional - An optional callback function (or callable) that is called with the + :type iterable: iterable + :param callback: An optional callback function (or callable) that is called with the result from each worker run and is executed on the master process. This is useful for, e.g., saving results to a file, since the callback is only called on the master thread. + :type callback: callable, optional - Returns - ------- - results : list - A list of results from the output of each ``worker()`` call. - + :return: A list of results from the output of each ``worker()`` call. """ if callback is None: diff --git a/lenstronomy/Sampling/Pool/pool.py b/lenstronomy/Sampling/Pool/pool.py index 837b9ea15..59256ed7b 100644 --- a/lenstronomy/Sampling/Pool/pool.py +++ b/lenstronomy/Sampling/Pool/pool.py @@ -50,18 +50,16 @@ def choose_pool(mpi=False, processes=1, **kwargs): Choose between the different pools given options from, e.g., argparse. - Parameters - ---------- - mpi : bool, optional - Use the MPI processing pool, :class:`~schwimmbad.mpi.MPIPool`. By + + :param mpi: Use the MPI processing pool, :class:`~schwimmbad.mpi.MPIPool`. By default, ``False``, will use the :class:`~schwimmbad.serial.SerialPool`. - processes : int, optional - Use the multiprocessing pool, + :type mpi: bool, optional + :param processes: Use the multiprocessing pool, :class:`~schwimmbad.multiprocessing.MultiPool`, with this number of processes. By default, ``processes=1``, will use them:class:`~schwimmbad.serial.SerialPool`. - - Any additional kwargs are passed in to the pool class initializer selected by the arguments. - + :type processes: int, optional + :param kwargs: Any additional kwargs are passed in to the pool class initializer selected by the arguments. + :type kwargs: keyword arguments """ # Imports moved here to avoid crashing at import time if dependencies # are missing diff --git a/lenstronomy/Sampling/Samplers/nautilus.py b/lenstronomy/Sampling/Samplers/nautilus.py index 70c9564a8..e6948827d 100644 --- a/lenstronomy/Sampling/Samplers/nautilus.py +++ b/lenstronomy/Sampling/Samplers/nautilus.py @@ -44,12 +44,11 @@ def nautilus_sampling(self, prior_type='uniform', mpi=False, thread_count=1, ver if prior_type == 'uniform': for i in range(self._num_param): prior.add_parameter(dist=(self._lower_limit[i], self._upper_limit[i])) - print(self._num_param, prior, 'test') else: raise ValueError('prior_type %s is not supported for Nautilus wrapper.' % prior_type) # loop through prior pool = choose_pool(mpi=mpi, processes=thread_count, use_dill=True) - sampler = Sampler(prior, likelihood=self.likelihood, pool=pool, **kwargs_nautilus) + sampler = Sampler(prior, likelihood=self.likelihood, pool=pool, pass_struct=False, **kwargs_nautilus) time_start = time.time() if one_step is True: sampler.add_bound() diff --git a/lenstronomy/Sampling/Samplers/pso.py b/lenstronomy/Sampling/Samplers/pso.py index e180cf5cf..f541e0aa3 100644 --- a/lenstronomy/Sampling/Samplers/pso.py +++ b/lenstronomy/Sampling/Samplers/pso.py @@ -105,6 +105,7 @@ def set_global_best(self, position, velocity, fitness): def _init_swarm(self): """ Initiate the swarm. + :return: :rtype: """ @@ -125,10 +126,9 @@ def sample(self, max_iter=1000, c1=1.193, c2=1.193, p=0.7, m=1e-3, n=1e-2, early :param c1: cognitive weight :param c2: social weight :param p: stop criterion, percentage of particles to use - :param m: stop criterion, difference between mean fitness and global - best + :param m: stop criterion, difference between mean fitness and global best :param n: stop criterion, difference between norm of the particle - vector and norm of the global best + vector and norm of the global best :param early_stop_tolerance: will terminate at the given value (should be specified as a chi^2) """ @@ -194,19 +194,18 @@ def optimize(self, max_iter=1000, verbose=True, c1=1.193, c2=1.193, :param c1: cognitive weight :param c2: social weight :param p: stop criterion, percentage of particles to use - :param m: stop criterion, difference between mean fitness and global - best + :param m: stop criterion, difference between mean fitness and global best :param n: stop criterion, difference between norm of the particle - vector and norm of the global best + vector and norm of the global best :param early_stop_tolerance: will terminate at the given value (should be specified as a chi^2) """ - chi2_list = [] + log_likelihood_list = [] vel_list = [] pos_list = [] num_iter = 0 for _ in self.sample(max_iter, c1, c2, p, m, n, early_stop_tolerance): - chi2_list.append(self.global_best.fitness * 2) + log_likelihood_list.append(self.global_best.fitness) vel_list.append(self.global_best.velocity) pos_list.append(self.global_best.position) num_iter += 1 @@ -215,11 +214,12 @@ def optimize(self, max_iter=1000, verbose=True, c1=1.193, c2=1.193, if num_iter % 10 == 0: print(num_iter) - return self.global_best.position, [chi2_list, pos_list, vel_list] + return self.global_best.position, [log_likelihood_list, pos_list, vel_list] def _get_fitness(self, swarm): """ Set fitness (probability) of the particles in swarm. + :param swarm: PSO state :type swarm: list of Particle() instances of the swarm :return: @@ -239,6 +239,7 @@ def _get_fitness(self, swarm): def _converged(self, it, p, m, n): """ Check for convergence. + :param it: :type it: :param p: @@ -274,7 +275,7 @@ def _converged_fit(self, it, p, m): best_sort = np.sort([particle.personal_best.fitness for particle in self.swarm])[::-1] mean_fit = np.mean(best_sort[1:int(math.floor(self.particleCount * p))]) - #print( "best %f, mean_fit %f, ration %f"%( self.global_best[0], + # print( "best %f, mean_fit %f, ration %f"%( self.global_best[0], # mean_fit, abs((self.global_best[0]-mean_fit)))) return abs(self.global_best.fitness - mean_fit) < m diff --git a/lenstronomy/Sampling/likelihood.py b/lenstronomy/Sampling/likelihood.py index 0b6b29fa4..f807c3cf0 100644 --- a/lenstronomy/Sampling/likelihood.py +++ b/lenstronomy/Sampling/likelihood.py @@ -43,7 +43,7 @@ def __init__(self, kwargs_data_joint, kwargs_model, param_class, image_likelihoo into the conventions of the imSim_class :param image_likelihood: bool, option to compute the imaging likelihood :param source_position_likelihood: bool, if True, ray-traces image positions back to source plane and evaluates - relative errors in respect ot the position_uncertainties in the image plane + relative errors in respect ot the position_uncertainties in the image plane :param check_bounds: bool, option to punish the hard bounds in parameter space :param check_matched_source_position: bool, option to check whether point source position of solver finds a solution to match all the image positions in the same source plane coordinate @@ -58,13 +58,13 @@ def __init__(self, kwargs_data_joint, kwargs_model, param_class, image_likelihoo :param image_likelihood_mask_list: list of boolean 2d arrays of size of images marking the pixels to be evaluated in the likelihood :param force_no_add_image: bool, if True: computes ALL image positions of the point source. If there are more - images predicted than modelled, a punishment occures + images predicted than modelled, a punishment occurs :param source_marg: marginalization addition on the imaging likelihood based on the covariance of the inferred linear coefficients :param linear_prior: float or list of floats (when multi-linear setting is chosen) indicating the range of - linear amplitude priors when computing the marginalization term. + linear amplitude priors when computing the marginalization term. :param restrict_image_number: bool, if True: computes ALL image positions of the point source. If there are more - images predicted than indicated in max_num_images, a punishment occurs + images predicted than indicated in max_num_images, a punishment occurs :param max_num_images: int, see restrict_image_number :param bands_compute: list of bools with same length as data objects, indicates which "band" to include in the fitting @@ -156,6 +156,13 @@ def __call__(self, a): def logL(self, args, verbose=False): """ routine to compute X2 given variable parameters for a MCMC/PSO chain + + + :param args: ordered parameter values that are being sampled + :type args: tuple or list of floats + :param verbose: if True, makes print statements about individual likelihood components + :type verbose: boolean + :returns: log likelihood of the data given the model (natural logarithm) """ # extract parameters kwargs_return = self.param.args2kwargs(args) @@ -166,6 +173,19 @@ def logL(self, args, verbose=False): return self.log_likelihood(kwargs_return, verbose=verbose) def log_likelihood(self, kwargs_return, verbose=False): + """ + + + :param kwargs_return: need to contain 'kwargs_lens', 'kwargs_source', 'kwargs_lens_light', 'kwargs_ps', + 'kwargs_special'. These entries themselves are lists of keyword argument of the parameters entering the model + to be evaluated + :type kwargs_return: keyword arguments + :param verbose: if True, makes print statements about individual likelihood components + :type verbose: boolean + + :returns: + - logL (float) log likelihood of the data given the model (natural logarithm) + """ kwargs_lens, kwargs_source, kwargs_lens_light, kwargs_ps, kwargs_special = kwargs_return['kwargs_lens'], \ kwargs_return['kwargs_source'], \ kwargs_return['kwargs_lens_light'], \ diff --git a/lenstronomy/Sampling/param_group.py b/lenstronomy/Sampling/param_group.py new file mode 100644 index 000000000..d85628622 --- /dev/null +++ b/lenstronomy/Sampling/param_group.py @@ -0,0 +1,340 @@ +''' +This module provides helper classes for managing sample parameters. This is +for internal use, if you are not modifying lenstronomy sampling to include +new parameters you can safely ignore this. +''' + +__author__ = 'jhodonnell' +__all__ = ['ModelParamGroup', 'SingleParam', 'ArrayParam'] + + +class ModelParamGroup: + ''' + This abstract class represents any lenstronomy fitting parameters used + in the Param class. + + Subclasses should implement num_params(), set_params(), and get_params() + to convert parameters from lenstronomy's semantic dictionary format to a + flattened array format and back. + + This class also contains three static methods to easily aggregate groups + of parameter classes, called `compose_num_params()`, `compose_set_params()`, + and `compose_get_params()`. + ''' + def num_params(self): + ''' + Tells the number of parameters that this group samples and their names. + + :returns: 2-tuple of (num param, list of names) + ''' + raise NotImplementedError + + def set_params(self, kwargs): + ''' + Converts lenstronomy semantic parameters in dictionary format into a + flattened array of parameters. + + The flattened array is for use in optimization algorithms, e.g. MCMC, + Particle swarm, etc. + + :returns: flattened array of parameters as floats + ''' + raise NotImplementedError + + def get_params(self, args, i): + ''' + Converts a flattened array of parameters back into a lenstronomy dictionary, + starting at index i. + + :param args: flattened arguments to convert to lenstronomy format + :type args: list + :param i: index to begin at in args + :type i: int + :returns: dictionary of parameters + ''' + raise NotImplementedError + + @staticmethod + def compose_num_params(each_group, *args, **kwargs): + ''' + Aggregates the number of parameters for a group of parameter groups, + calling each instance's `num_params()` method and combining the results + + :param each_group: collection of parameter groups. Should each be subclasses of ModelParamGroup. + :type each_group: list + :param args: Extra arguments to be passed to each call of `num_params()` + :param kwargs: Extra keyword arguments to be passed to each call of `num_params()` + + :returns: As in each individual `num_params()`, a 2-tuple of (num params, list of param names) + ''' + tot_param = 0 + param_names = [] + for group in each_group: + npar, names = group.num_params(*args, **kwargs) + tot_param += npar + param_names += names + return tot_param, param_names + + @staticmethod + def compose_set_params(each_group, param_kwargs, *args, **kwargs): + ''' + Converts lenstronomy semantic arguments in dictionary format to a + flattened list of floats for use in optimization/fitting algorithms. + Combines the results for a set of arbitrarily many parameter groups. + + :param each_group: collection of parameter groups. Should each be subclasses of ModelParamGroup. + :type each_group: list + :param param_kwargs: the kwargs to process + :type param_kwargs: dict + :param args: Extra arguments to be passed to each call of `set_params()` + :param kwargs: Extra keyword arguments to be passed to each call of `set_params()` + + :returns: As in each individual `set_params()`, a list of floats + ''' + output_args = [] + for group in each_group: + output_args += group.set_params(param_kwargs, *args, **kwargs) + return output_args + + @staticmethod + def compose_get_params(each_group, flat_args, i, *args, **kwargs): + ''' + Converts a flattened array of parameters to lenstronomy semantic + parameters in dictionary format. + Combines the results for a set of arbitrarily many parameter groups. + + :param each_group: collection of parameter groups. Should each be subclasses of ModelParamGroup. + :type each_group: list + :param flat_args: the input array of parameters + :type flat_args: list + :param i: the index in `flat_args` to start at + :type i: int + :param args: Extra arguments to be passed to each call of `set_params()` + :param kwargs: Extra keyword arguments to be passed to each call of `set_params()` + + :returns: As in each individual `get_params()`, a 2-tuple of (dictionary of params, new index) + ''' + output_kwargs = {} + for group in each_group: + kwargs_grp, i = group.get_params(flat_args, i, *args, **kwargs) + output_kwargs = dict(output_kwargs, **kwargs_grp) + return output_kwargs, i + + +class SingleParam(ModelParamGroup): + ''' + Helper for handling parameters which are a single float. + + Subclasses should define: + + :param on: Whether this parameter is sampled + :type on: bool + :param param_names: List of strings, the name of each parameter + :param _kwargs_lower: Dictionary. Lower bounds of each parameter + :param _kwargs_upper: Dictionary. Upper bounds of each parameter + ''' + def __init__(self, on): + ''' + :param on: Whether this paramter should be sampled + :type on: bool + ''' + self._on = bool(on) + + def num_params(self, kwargs_fixed): + ''' + Tells the number of parameters that this group samples and their names. + + :param kwargs_fixed: Dictionary of fixed arguments + :type kwargs_fixed: dict + + :returns: 2-tuple of (num param, list of names) + ''' + if self.on: + npar, names = 0, [] + for name in self.param_names: + if name not in kwargs_fixed: + npar += 1 + names.append(name) + return npar, names + return 0, [] + + def set_params(self, kwargs, kwargs_fixed): + ''' + Converts lenstronomy semantic parameters in dictionary format into a + flattened array of parameters. + + The flattened array is for use in optimization algorithms, e.g. MCMC, + Particle swarm, etc. + + :param kwargs: lenstronomy parameters to flatten + :type kwargs: dict + :param kwargs_fixed: Dictionary of fixed arguments + :type kwargs_fixed: dict + + :returns: flattened array of parameters as floats + ''' + if self.on: + output = [] + for name in self.param_names: + if name not in kwargs_fixed: + output.append(kwargs[name]) + return output + return [] + + def get_params(self, args, i, kwargs_fixed): + ''' + Converts a flattened array of parameters back into a lenstronomy dictionary, + starting at index i. + + :param args: flattened arguments to convert to lenstronomy format + :type args: list + :param i: index to begin at in args + :type i: int + :param kwargs_fixed: Dictionary of fixed arguments + :type kwargs_fixed: dict + + :returns: dictionary of parameters + ''' + out = {} + if self.on: + for name in self.param_names: + if name in kwargs_fixed: + out[name] = kwargs_fixed[name] + else: + out[name] = args[i] + i += 1 + return out, i + + @property + def kwargs_lower(self): + if not self.on: + return {} + return self._kwargs_lower + + @property + def kwargs_upper(self): + if not self.on: + return {} + return self._kwargs_upper + + @property + def on(self): + return self._on + + +class ArrayParam(ModelParamGroup): + ''' + Helper for handling parameters which are an array of values. Examples + include mass_scaling, which is an array of scaling parameters, and wavelet + or gaussian decompositions which have different coefficients for each mode. + + Subclasses should define: + + :param on: Whether this parameter is sampled + :type on: bool + :param param_names: Dictionary mapping the name of each parameter to the number of values needed. + :param _kwargs_lower: Dictionary. Lower bounds of each parameter + :param _kwargs_upper: Dictionary. Upper bounds of each parameter + ''' + def __init__(self, on): + ''' + :param on: Whether this paramter should be sampled + :type on: bool + ''' + self._on = bool(on) + + def num_params(self, kwargs_fixed): + ''' + Tells the number of parameters that this group samples and their names. + + :param kwargs_fixed: Dictionary of fixed arguments + :type kwargs_fixed: dict + + :returns: 2-tuple of (num param, list of names) + ''' + if not self.on: + return 0, [] + + npar = 0 + names = [] + for name, count in self.param_names.items(): + if name not in kwargs_fixed: + npar += count + names += [name] * count + + return npar, names + + def set_params(self, kwargs, kwargs_fixed): + ''' + Converts lenstronomy semantic parameters in dictionary format into a + flattened array of parameters. + + The flattened array is for use in optimization algorithms, e.g. MCMC, + Particle swarm, etc. + + :param kwargs: lenstronomy parameters to flatten + :type kwargs: dict + :param kwargs_fixed: Dictionary of fixed arguments + :type kwargs_fixed: dict + + :returns: flattened array of parameters as floats + ''' + if not self.on: + return [] + + args = [] + for name, count in self.param_names.items(): + if name not in kwargs_fixed: + args.extend(kwargs[name]) + return args + + def get_params(self, args, i, kwargs_fixed): + ''' + Converts a flattened array of parameters back into a lenstronomy dictionary, + starting at index i. + + :param args: flattened arguments to convert to lenstronomy format + :type args: list + :param i: index to begin at in args + :type i: int + :param kwargs_fixed: Dictionary of fixed arguments + :type kwargs_fixed: dict + + :returns: dictionary of parameters + ''' + if not self.on: + return {}, i + + params = {} + for name, count in self.param_names.items(): + if name not in kwargs_fixed: + params[name] = args[i:i + count] + i += count + else: + params[name] = kwargs_fixed[name] + + return params, i + + @property + def kwargs_lower(self): + if not self.on: + return {} + + out = {} + for name, count in self.param_names.items(): + out[name] = [self._kwargs_lower[name]] * count + return out + + @property + def kwargs_upper(self): + if not self.on: + return {} + + out = {} + for name, count in self.param_names.items(): + out[name] = [self._kwargs_upper[name]] * count + return out + + @property + def on(self): + return self._on diff --git a/lenstronomy/Sampling/parameters.py b/lenstronomy/Sampling/parameters.py index 45640dc39..3cab1acc7 100644 --- a/lenstronomy/Sampling/parameters.py +++ b/lenstronomy/Sampling/parameters.py @@ -47,6 +47,48 @@ class that handles the parameter constraints. In particular when different model 'joint_lens_with_source_light': [[i_source, k_lens, ['param_name1', 'param_name2', ...]], [...], ...], joint parameter between lens model and source light model. Samples light model parameter only. + 'mass_scaling_list': e.g. [False, 1, 1, False, 2, False, 1, ...] + Links lens models to have their masses scaled together. In this example, + masses with False are not scaled, masses labeled 1 are scaled together, + and those labeled 2 are scaled together independent of 1, etc. + + 'general_scaling': { 'param1': [False, 1, 1, False, 1, ...], 'param2': [1, 1, 1, False, 2, 2, ...] } + Generalized parameter scaling. Input should be a dictionary mapping + parameter names to the masks defining which lens models are scaled together, + in the same format as for 'mass_scaling_list'. For each scaled parameter, + two special params will be added called '${param}_scale_factor' and + '${param}_scale_pow', defining the scaling and power-law of each. + + Each scale will be modified as `param = param_scale_factor * param**param_scale_pow`. + + For example, say we want to jointly constrain the `sigma0` and `Rs` parameters + of some lens models indexed by `i`, like so: + + .. math:: + + \\sigma_{0,i} = \\sigma_0^{ref} L_i^\\alpha \\\\ + r_{cut,i} = r_{cut}^{ref} L_i^\\beta + + To do this we can add the following. The lens models corresponding to + entries of `1` will be scaled together, and those corresponding to `False` + will not be. As in `mass_scaling_list`, subsets of models can be scaled + independently by marking them `2`, `3`, etc. + + >>> 'general_scaling': { + >>> 'sigma0': [False, 1, 1, False, 1, ...], + >>> 'Rs': [False, 1, 1, False, 1, ...], + >>> } + + Then we can choose to fix the power-law and vary the scale factor like so: + + >>> fixed_special = {'sigma0_scale_pow': [alpha*2], 'Rs_scale_pow': [beta]} + >>> kwargs_special_init = {'sigma0_scale_factor': [17.0], 'Rs_scale_factor': [8]} + >>> kwargs_special_sigma = {'sigma0_scale_factor': [10.0], 'Rs_scale_factor': [3]} + >>> kwargs_lower_special = {'sigma0_scale_factor': [0.5], 'Rs_scale_factor': [1]} + >>> kwargs_upper_special = {'sigma0_scale_factor': [40], 'Rs_scale_factor': [20]} + + >>> special_params = [kwargs_special_init, kwargs_special_sigma, fixed_special, kwargs_lower_special, kwargs_upper_special] + hierarchy is as follows: 1. Point source parameters are inferred 2. Lens light joint parameters are set @@ -77,6 +119,7 @@ def __init__(self, kwargs_model, joint_source_with_source=[], joint_lens_with_light=[], joint_source_with_point_source=[], joint_lens_light_with_point_source=[], joint_extinction_with_lens_light=[], joint_lens_with_source_light=[], mass_scaling_list=None, point_source_offset=False, + general_scaling=None, num_point_source_list=None, image_plane_source_list=None, solver_type='NONE', Ddt_sampling=None, source_size=False, num_tau0=0, lens_redshift_sampling_indexes=None, source_redshift_sampling_indexes=None, source_grid_offset=False, num_shapelet_lens=0, @@ -125,6 +168,7 @@ def __init__(self, kwargs_model, joint parameter between lens model and source light model. Samples light model parameter only. :param mass_scaling_list: boolean list of length of lens model list (optional) models with identical integers will be scaled with the same additional scaling factor. First integer starts with 1 (not 0) + :param general_scaling: { 'param_1': [list of booleans/integers defining which model to fit], 'param_2': [..], ..} :param point_source_offset: bool, if True, adds relative offsets ot the modeled image positions relative to the time-delay and lens equation solver :param num_point_source_list: list of number of point sources per point source model class @@ -144,7 +188,7 @@ def __init__(self, kwargs_model, in ascending numbering e.g. [-1, 0, 0, 1, 0, 2], -1 indicating not sampled fixed indexes. These indexes are the sample as for the lens :param source_grid_offset: optional, if True when using a pixel-based modelling (e.g. with STARLETS-like profiles), - adds two additional sampled parameters describing RA/Dec offsets between data coordinate grid and pixelated source plane coordinate grid. + adds two additional sampled parameters describing RA/Dec offsets between data coordinate grid and pixelated source plane coordinate grid. :param num_shapelet_lens: number of shapelet coefficients in the 'SHAPELETS_CART' or 'SHAPELETS_POLAR' mass profile. :param log_sampling_lens: Sample the log10 of the lens model parameters. Format : [[i_lens, ['param_name1', 'param_name2', ...]], [...], ...], """ @@ -211,6 +255,18 @@ def __init__(self, kwargs_model, else: self._num_scale_factor = 0 self._mass_scaling = False + + if general_scaling is not None: + self._general_scaling = True + # FIXME TODO: check that the scaled parameters are actually used + # by each lens model. For example, NFW has no theta_E parameter, + # if a user tries to scale theta_E for an NFW we should throw an + # error + self._general_scaling_masks = dict(general_scaling) + else: + self._general_scaling = False + self._general_scaling_masks = dict() + self._point_source_offset = point_source_offset if num_point_source_list is None: num_point_source_list = [1] * len(self._point_source_model_list) @@ -270,6 +326,7 @@ def __init__(self, kwargs_model, kwargs_lower=kwargs_lower_extinction, kwargs_upper=kwargs_upper_extinction, linear_solver=False) self.specialParams = SpecialParam(Ddt_sampling=Ddt_sampling, mass_scaling=self._mass_scaling, + general_scaling_params=self._general_scaling_masks, kwargs_fixed=kwargs_fixed_special, num_scale_factor=self._num_scale_factor, kwargs_lower=kwargs_lower_special, kwargs_upper=kwargs_upper_special, point_source_offset=self._point_source_offset, num_images=self._num_images, @@ -544,24 +601,43 @@ def update_lens_scaling(self, kwargs_special, kwargs_lens, inverse=False): :return: updated lens model keyword argument list """ kwargs_lens_updated = copy.deepcopy(kwargs_lens) - if self._mass_scaling is False: + # If we do not scaling, there's nothing to be done + if not (self._mass_scaling or self._general_scaling): return kwargs_lens_updated - scale_factor_list = np.array(kwargs_special['scale_factor']) - if inverse is True: - scale_factor_list = 1. / np.array(kwargs_special['scale_factor']) - for i, kwargs in enumerate(kwargs_lens_updated): - if self._mass_scaling_list[i] is not False: - scale_factor = scale_factor_list[self._mass_scaling_list[i] - 1] - if 'theta_E' in kwargs: - kwargs['theta_E'] *= scale_factor - elif 'alpha_Rs' in kwargs: - kwargs['alpha_Rs'] *= scale_factor - elif 'alpha_1' in kwargs: - kwargs['alpha_1'] *= scale_factor - elif 'sigma0' in kwargs: - kwargs['sigma0'] *= scale_factor - elif 'k_eff' in kwargs: - kwargs['k_eff'] *= scale_factor + + # TODO: remove separate logic for mass scaling. either deprecate it + # entirely, implement the details as a special case of general_scaling + if self._mass_scaling: + scale_factor_list = np.array(kwargs_special['scale_factor']) + if inverse is True: + scale_factor_list = 1. / np.array(kwargs_special['scale_factor']) + for i, kwargs in enumerate(kwargs_lens_updated): + if self._mass_scaling_list[i] is not False: + scale_factor = scale_factor_list[self._mass_scaling_list[i] - 1] + if 'theta_E' in kwargs: + kwargs['theta_E'] *= scale_factor + elif 'alpha_Rs' in kwargs: + kwargs['alpha_Rs'] *= scale_factor + elif 'alpha_1' in kwargs: + kwargs['alpha_1'] *= scale_factor + elif 'sigma0' in kwargs: + kwargs['sigma0'] *= scale_factor + elif 'k_eff' in kwargs: + kwargs['k_eff'] *= scale_factor + + if self._general_scaling: + for param_name in self._general_scaling_masks.keys(): + factors = kwargs_special[f'{param_name}_scale_factor'] + _pows = kwargs_special[f'{param_name}_scale_pow'] + + for i, kwargs in enumerate(kwargs_lens_updated): + scale_idx = self._general_scaling_masks[param_name][i] + if scale_idx is not False: + if inverse: + kwargs[param_name] = (kwargs[param_name] / factors[scale_idx - 1]) ** (1 / _pows[scale_idx - 1]) + else: + kwargs[param_name] = factors[scale_idx - 1] * kwargs[param_name]**_pows[scale_idx - 1] + return kwargs_lens_updated def _add_fixed_lens(self, kwargs_fixed, kwargs_init): @@ -642,6 +718,7 @@ def print_setting(self): num, param_list = self.num_param() num_linear = self.num_param_linear() + # TODO print settings of specailParams? print("The following model options are chosen:") print("Lens models:", self._lens_model_list) print("Source models:", self._source_light_model_list) @@ -661,6 +738,8 @@ def print_setting(self): print("Joint lens with light:", self._joint_lens_with_light) print("Joint source with point source:", self._joint_source_with_point_source) print("Joint lens light with point source:", self._joint_lens_light_with_point_source) + print("Mass scaling:", self._num_scale_factor, "groups") + print("General lens scaling:", self._general_scaling_masks) print("===================") print("Number of non-linear parameters being sampled: ", num) print("Parameters being sampled: ", param_list) diff --git a/lenstronomy/Sampling/sampler.py b/lenstronomy/Sampling/sampler.py index 9e68b2e92..7e0ffd672 100644 --- a/lenstronomy/Sampling/sampler.py +++ b/lenstronomy/Sampling/sampler.py @@ -44,7 +44,7 @@ def simplex(self, init_pos, n_iterations, method, print_key='SIMPLEX'): kwargs_return = self.chain.param.args2kwargs(result['x']) print(-logL * 2 / (max(self.chain.effective_num_data_points(**kwargs_return), 1)), 'reduced X^2 of best position') - print(logL, 'logL') + print(logL, 'log likelihood') print(self.chain.effective_num_data_points(**kwargs_return), 'effective number of data points') print(kwargs_return.get('kwargs_lens', None), 'lens result') print(kwargs_return.get('kwargs_source', None), 'source result') @@ -100,13 +100,13 @@ def pso(self, n_particles, n_iterations, lower_start=None, upper_start=None, time_start = time.time() - result, [chi2_list, pos_list, vel_list] = pso.optimize(n_iterations) + result, [log_likelihood_list, pos_list, vel_list] = pso.optimize(n_iterations) if pool.is_master(): kwargs_return = self.chain.param.args2kwargs(result) print(pso.global_best.fitness * 2 / (max( self.chain.effective_num_data_points(**kwargs_return), 1)), 'reduced X^2 of best position') - print(pso.global_best.fitness, 'logL') + print(pso.global_best.fitness, 'log likelihood') print(self.chain.effective_num_data_points(**kwargs_return), 'effective number of data points') print(kwargs_return.get('kwargs_lens', None), 'lens result') print(kwargs_return.get('kwargs_source', None), 'source result') @@ -116,7 +116,7 @@ def pso(self, n_particles, n_iterations, lower_start=None, upper_start=None, time_end = time.time() print(time_end - time_start, 'time used for ', print_key) print('===================') - return result, [chi2_list, pos_list, vel_list] + return result, [log_likelihood_list, pos_list, vel_list] def mcmc_emcee(self, n_walkers, n_run, n_burn, mean_start, sigma_start, mpi=False, progress=False, threadCount=1, @@ -217,12 +217,14 @@ def mcmc_zeus(self, n_walkers, n_run, n_burn, mean_start, sigma_start, :type mean_start: numpy array of length the number of parameters :param sigma_start: spread of the parameter values (uncorrelated in each dimension) of the initialising sample :type sigma_start: numpy array of length the number of parameters + :param mpi: if True, initializes an MPIPool to allow for MPI execution of the sampler + :type mpi: bool :param progress: :type progress: bool :param initpos: initial walker position to start sampling (optional) :type initpos: numpy array of size num param x num walkser - :param backup_filename: name of the HDF5 file where sampling state is saved (through zeus callback function) - :type backup_filename: string + :param backend_filename: name of the HDF5 file where sampling state is saved (through zeus callback function) + :type backend_filename: string :return: samples, ln likelihood value of samples :rtype: numpy 2d array, numpy 1d array """ @@ -251,7 +253,7 @@ def mcmc_zeus(self, n_walkers, n_run, n_burn, mean_start, sigma_start, blobs_dtype=blobs_dtype, verbose=verbose, check_walkers=check_walkers, shuffle_ensemble=shuffle_ensemble, light_mode=light_mode) - sampler.run_mcmc(initpos, n_run_eff, progress=progress, callbacks = backend) + sampler.run_mcmc(initpos, n_run_eff, progress=progress, callbacks=backend) flat_samples = sampler.get_chain(flat=True, thin=1, discard=n_burn) diff --git a/lenstronomy/Sampling/special_param.py b/lenstronomy/Sampling/special_param.py index 95d5d355e..7cc8e7391 100644 --- a/lenstronomy/Sampling/special_param.py +++ b/lenstronomy/Sampling/special_param.py @@ -2,6 +2,133 @@ __all__ = ['SpecialParam'] +import numpy as np +from .param_group import ModelParamGroup, SingleParam, ArrayParam + + +# ==================================== # +# == Defining individual parameters == # +# ==================================== # + + +class DdtSamplingParam(SingleParam): + ''' + Time delay parameter + ''' + param_names = ['D_dt'] + _kwargs_lower = {'D_dt': 0} + _kwargs_upper = {'D_dt': 100000} + + +class SourceSizeParam(SingleParam): + ''' + Source size parameter + ''' + param_names = ['source_size'] + _kwargs_lower = {'source_size': 0} + _kwargs_upper = {'source_size': 1} + + +class SourceGridOffsetParam(SingleParam): + ''' + Source grid offset, both x and y. + ''' + param_names = ['delta_x_source_grid', 'delta_y_source_grid'] + _kwargs_lower = { + 'delta_x_source_grid': -100, + 'delta_y_source_grid': -100 + } + _kwargs_upper = { + 'delta_x_source_grid': 100, + 'delta_y_source_grid': 100 + } + + +class MassScalingParam(ArrayParam): + ''' + Mass scaling. Can scale the masses of arbitrary subsets of lens models + ''' + _kwargs_lower = {'scale_factor': 0} + _kwargs_upper = {'scale_factor': 1000} + def __init__(self, num_scale_factor): + super().__init__(on=int(num_scale_factor) > 0) + self.param_names = {'scale_factor': int(num_scale_factor)} + + +class PointSourceOffsetParam(ArrayParam): + ''' + Point source offset, both x and y + ''' + _kwargs_lower = {'delta_x_image': -1, 'delta_y_image': -1} + _kwargs_upper = {'delta_x_image': 1, 'delta_y_image': 1} + def __init__(self, offset, num_images): + super().__init__(on=offset and (int(num_images) > 0)) + self.param_names = { + 'delta_x_image': int(num_images), + 'delta_y_image': int(num_images), + } + + +class Tau0ListParam(ArrayParam): + ''' + Optical depth renormalization parameters + ''' + _kwargs_lower = {'tau0_list': 0} + _kwargs_upper = {'tau0_list': 1000} + def __init__(self, num_tau0): + super().__init__(on=int(num_tau0) > 0) + self.param_names = {'tau0_list': int(num_tau0)} + + +class ZSamplingParam(ArrayParam): + ''' + Redshift sampling. + ''' + _kwargs_lower = {'z_sampling': 0} + _kwargs_upper = {'z_sampling': 1000} + def __init__(self, num_z_sampling): + super().__init__(on=int(num_z_sampling) > 0) + self.param_names = {'z_sampling': int(num_z_sampling)} + + +class GeneralScalingParam(ArrayParam): + ''' + General lens scaling. + + For each scaled lens parameter, adds a `{param}_scale_factor` and + `{param}_scale_pow` special parameter, and updates the scaled param + as `param = param_scale_factor * param**param_scale_pow`. + ''' + def __init__(self, params: dict): + # params is a dictionary + self.param_names = {} + self._kwargs_lower = {} + self._kwargs_upper = {} + + super().__init__(params) + if not self.on: + return + + for name, array in params.items(): + num_param = np.max(array) + + if num_param > 0: + fac_name = f'{name}_scale_factor' + self.param_names[fac_name] = num_param + self._kwargs_lower[fac_name] = 0 + self._kwargs_upper[fac_name] = 1000 + + pow_name = f'{name}_scale_pow' + self.param_names[pow_name] = num_param + self._kwargs_lower[pow_name] = -10 + self._kwargs_upper[pow_name] = 10 + + +# ======================================== # +# == All together: Composing into class == # +# ======================================== # + + class SpecialParam(object): """ @@ -9,7 +136,8 @@ class that handles special parameters that are not directly part of a specific m These includes cosmology relevant parameters, astrometric errors and overall scaling parameters. """ - def __init__(self, Ddt_sampling=False, mass_scaling=False, num_scale_factor=1, kwargs_fixed=None, kwargs_lower=None, + def __init__(self, Ddt_sampling=False, mass_scaling=False, num_scale_factor=1, + general_scaling_params=None, kwargs_fixed=None, kwargs_lower=None, kwargs_upper=None, point_source_offset=False, source_size=False, num_images=0, num_tau0=0, num_z_sampling=0, source_grid_offset=False): """ @@ -32,59 +160,35 @@ def __init__(self, Ddt_sampling=False, mass_scaling=False, num_scale_factor=1, k Warning: this is only defined for pixel-based source modelling (e.g. 'SLIT_STARLETS' light profile) """ - self._D_dt_sampling = Ddt_sampling - self._mass_scaling = mass_scaling - self._num_scale_factor = num_scale_factor - self._point_source_offset = point_source_offset - self._num_images = num_images - self._num_tau0 = num_tau0 - self._num_z_sampling = num_z_sampling - if num_z_sampling > 0: - self._z_sampling = True + self._D_dt_sampling = DdtSamplingParam(Ddt_sampling) + if not mass_scaling: + num_scale_factor = 0 + self._mass_scaling = MassScalingParam(num_scale_factor) + + self._general_scaling = GeneralScalingParam(general_scaling_params or dict()) + + if point_source_offset: + self._point_source_offset = PointSourceOffsetParam(True, num_images) else: - self._z_sampling = False + self._point_source_offset = PointSourceOffsetParam(False, 0) + self._source_size = SourceSizeParam(source_size) + self._tau0 = Tau0ListParam(num_tau0) + self._z_sampling = ZSamplingParam(num_z_sampling) + self._source_grid_offset = SourceGridOffsetParam(source_grid_offset) if kwargs_fixed is None: kwargs_fixed = {} self._kwargs_fixed = kwargs_fixed - self._source_size = source_size - self._source_grid_offset = source_grid_offset + if kwargs_lower is None: kwargs_lower = {} - if self._D_dt_sampling is True: - kwargs_lower['D_dt'] = 0 - if self._mass_scaling is True: - kwargs_lower['scale_factor'] = [0] * self._num_scale_factor - if self._point_source_offset is True: - kwargs_lower['delta_x_image'] = [-1] * self._num_images - kwargs_lower['delta_y_image'] = [-1] * self._num_images - if self._source_size is True: - kwargs_lower['source_size'] = 0 - if self._num_tau0 > 0: - kwargs_lower['tau0_list'] = [0] * self._num_tau0 - if self._z_sampling is True: - kwargs_lower['z_sampling'] = [0] * self._num_z_sampling - if self._source_grid_offset: - kwargs_lower['delta_x_source_grid'] = -100 - kwargs_lower['delta_y_source_grid'] = -100 + for group in self._param_groups: + kwargs_lower = dict(kwargs_lower, **group.kwargs_lower) if kwargs_upper is None: kwargs_upper = {} - if self._D_dt_sampling is True: - kwargs_upper['D_dt'] = 100000 - if self._mass_scaling is True: - kwargs_upper['scale_factor'] = [1000] * self._num_scale_factor - if self._point_source_offset is True: - kwargs_upper['delta_x_image'] = [1] * self._num_images - kwargs_upper['delta_y_image'] = [1] * self._num_images - if self._source_size is True: - kwargs_upper[source_size] = 1 - if self._num_tau0 > 0: - kwargs_upper['tau0_list'] = [1000] * self._num_tau0 - if self._z_sampling is True: - kwargs_upper['z_sampling'] = [20] * self._num_z_sampling - if self._source_grid_offset: - kwargs_upper['delta_x_source_grid'] = 100 - kwargs_upper['delta_y_source_grid'] = 100 + for group in self._param_groups: + kwargs_upper = dict(kwargs_upper, **group.kwargs_upper) + self.lower_limit = kwargs_lower self.upper_limit = kwargs_upper @@ -95,60 +199,10 @@ def get_params(self, args, i): :param i: integer, list index to start the read out for this class :return: keyword arguments related to args, index after reading out arguments of this class """ - kwargs_special = {} - if self._D_dt_sampling is True: - if 'D_dt' not in self._kwargs_fixed: - kwargs_special['D_dt'] = args[i] - i += 1 - else: - kwargs_special['D_dt'] = self._kwargs_fixed['D_dt'] - if self._mass_scaling is True: - if 'scale_factor' not in self._kwargs_fixed: - kwargs_special['scale_factor'] = args[i: i + self._num_scale_factor] - i += self._num_scale_factor - else: - kwargs_special['scale_factor'] = self._kwargs_fixed['scale_factor'] - if self._point_source_offset is True: - if 'delta_x_image' not in self._kwargs_fixed: - kwargs_special['delta_x_image'] = args[i: i + self._num_images] - i += self._num_images - else: - kwargs_special['delta_x_image'] = self._kwargs_fixed['delta_x_image'] - if 'delta_y_image' not in self._kwargs_fixed: - kwargs_special['delta_y_image'] = args[i: i + self._num_images] - i += self._num_images - else: - kwargs_special['delta_y_image'] = self._kwargs_fixed['delta_y_image'] - if self._source_size is True: - if 'source_size' not in self._kwargs_fixed: - kwargs_special['source_size'] = args[i] - i += 1 - else: - kwargs_special['source_size'] = self._kwargs_fixed['source_size'] - if self._num_tau0 > 0: - if 'tau0_list' not in self._kwargs_fixed: - kwargs_special['tau0_list'] = args[i:i + self._num_tau0] - i += self._num_tau0 - else: - kwargs_special['tau0_list'] = self._kwargs_fixed['tau0_list'] - if self._z_sampling is True: - if 'z_sampling' not in self._kwargs_fixed: - kwargs_special['z_sampling'] = args[i:i + self._num_z_sampling] - i += self._num_z_sampling - else: - kwargs_special['z_sampling'] = self._kwargs_fixed['z_sampling'] - if self._source_grid_offset: - if 'delta_x_source_grid' not in self._kwargs_fixed: - kwargs_special['delta_x_source_grid'] = args[i] - i += 1 - else: - kwargs_special['delta_x_source_grid'] = self._kwargs_fixed['delta_x_source_grid'] - if 'delta_y_source_grid' not in self._kwargs_fixed: - kwargs_special['delta_y_source_grid'] = args[i] - i += 1 - else: - kwargs_special['delta_y_source_grid'] = self._kwargs_fixed['delta_y_source_grid'] - return kwargs_special, i + result = ModelParamGroup.compose_get_params( + self._param_groups, args, i, kwargs_fixed=self._kwargs_fixed + ) + return result def set_params(self, kwargs_special): """ @@ -156,83 +210,26 @@ def set_params(self, kwargs_special): :param kwargs_special: keyword arguments with parameter settings :return: argument list of the sampled parameters extracted from kwargs_special """ - args = [] - if self._D_dt_sampling is True: - if 'D_dt' not in self._kwargs_fixed: - args.append(kwargs_special['D_dt']) - if self._mass_scaling is True: - if 'scale_factor' not in self._kwargs_fixed: - for i in range(self._num_scale_factor): - args.append(kwargs_special['scale_factor'][i]) - if self._point_source_offset is True: - if 'delta_x_image' not in self._kwargs_fixed: - for i in range(self._num_images): - args.append(kwargs_special['delta_x_image'][i]) - if 'delta_y_image' not in self._kwargs_fixed: - for i in range(self._num_images): - args.append(kwargs_special['delta_y_image'][i]) - if self._source_size is True: - if 'source_size' not in self._kwargs_fixed: - args.append(kwargs_special['source_size']) - if self._num_tau0 > 0: - if 'tau0_list' not in self._kwargs_fixed: - for i in range(self._num_tau0): - args.append(kwargs_special['tau0_list'][i]) - if self._z_sampling is True: - if 'z_sampling' not in self._kwargs_fixed: - for i in range(self._num_z_sampling): - args.append(kwargs_special['z_sampling'][i]) - if self._source_grid_offset is True: - if 'delta_x_source_grid' not in self._kwargs_fixed: - args.append(kwargs_special['delta_x_source_grid']) - if 'delta_y_source_grid' not in self._kwargs_fixed: - args.append(kwargs_special['delta_y_source_grid']) - return args + return ModelParamGroup.compose_set_params( + self._param_groups, kwargs_special, kwargs_fixed=self._kwargs_fixed + ) def num_param(self): """ :return: integer, number of free parameters sampled (and managed) by this class, parameter names (list of strings) """ - num = 0 - string_list = [] - if self._D_dt_sampling is True: - if 'D_dt' not in self._kwargs_fixed: - num += 1 - string_list.append('D_dt') - if self._mass_scaling is True: - if 'scale_factor' not in self._kwargs_fixed: - num += self._num_scale_factor - for i in range(self._num_scale_factor): - string_list.append('scale_factor') - if self._point_source_offset is True: - if 'delta_x_image' not in self._kwargs_fixed: - num += self._num_images - for i in range(self._num_images): - string_list.append('delta_x_image') - if 'delta_y_image' not in self._kwargs_fixed: - num += self._num_images - for i in range(self._num_images): - string_list.append('delta_y_image') - if self._source_size is True: - if 'source_size' not in self._kwargs_fixed: - num += 1 - string_list.append('source_size') - if self._num_tau0 > 0: - if 'tau0_list' not in self._kwargs_fixed: - num += self._num_tau0 - for i in range(self._num_tau0): - string_list.append('tau0') - if self._z_sampling is True: - if 'z_sampling' not in self._kwargs_fixed: - num += self._num_z_sampling - for i in range(self._num_z_sampling): - string_list.append('z') - if self._source_grid_offset is True: - if 'delta_x_source_grid' not in self._kwargs_fixed: - num += 1 - string_list.append('delta_x_source_grid') - if 'delta_y_source_grid' not in self._kwargs_fixed: - num += 1 - string_list.append('delta_y_source_grid') - return num, string_list + return ModelParamGroup.compose_num_params( + self._param_groups, kwargs_fixed=self._kwargs_fixed + ) + + @property + def _param_groups(self): + return [self._D_dt_sampling, + self._mass_scaling, + self._general_scaling, + self._point_source_offset, + self._source_size, + self._tau0, + self._z_sampling, + self._source_grid_offset] diff --git a/lenstronomy/SimulationAPI/data_api.py b/lenstronomy/SimulationAPI/data_api.py index 0e9a9d4a6..2fa0e148b 100644 --- a/lenstronomy/SimulationAPI/data_api.py +++ b/lenstronomy/SimulationAPI/data_api.py @@ -14,13 +14,21 @@ class DataAPI(SingleBand): options are available. Have a look in the specific modules if you are interested in. """ - def __init__(self, numpix, **kwargs_single_band): + def __init__(self, numpix, kwargs_pixel_grid=None, **kwargs_single_band): """ :param numpix: number of pixels per axis in the simulation to be modelled + :param kwargs_pixel_grid: if None, uses default pixel grid option + if defined, must contain keyword arguments PixelGrid() class :param kwargs_single_band: keyword arguments used to create instance of SingleBand class """ self.numpix = numpix + if kwargs_pixel_grid is not None: + required_keys = ['ra_at_xy_0','dec_at_xy_0','transform_pix2angle'] + if not all(k in kwargs_pixel_grid for k in required_keys): + raise ValueError('Missing 1 or more required'+ + 'kwargs_pixel_grid parameters') + self._kwargs_pixel_grid = kwargs_pixel_grid SingleBand.__init__(self, **kwargs_single_band) @property @@ -39,13 +47,22 @@ def kwargs_data(self): :return: keyword arguments for ImageData class instance """ - x_grid, y_grid, ra_at_xy_0, dec_at_xy_0, x_at_radec_0, y_at_radec_0, transform_pix2angle, transform_angle2pix = util.make_grid_with_coordtransform( - numPix=self.numpix, deltapix=self.pixel_scale, subgrid_res=1, left_lower=False, inverse=False) + # default pixel grid + if self._kwargs_pixel_grid is None: + _, _, ra_at_xy_0, dec_at_xy_0, _, _, transform_pix2angle, _ = util.make_grid_with_coordtransform( + numPix=self.numpix, deltapix=self.pixel_scale, subgrid_res=1, + left_lower=False, inverse=False) + # user defined pixel grid + else: + ra_at_xy_0 = self._kwargs_pixel_grid['ra_at_xy_0'] + dec_at_xy_0 = self._kwargs_pixel_grid['dec_at_xy_0'] + transform_pix2angle = self._kwargs_pixel_grid['transform_pix2angle'] # CCD gain corrected exposure time to allow a direct Poisson estimates based on IID counts scaled_exposure_time = self.flux_iid(1) - kwargs_data = {'image_data': np.zeros((self.numpix, self.numpix)), 'ra_at_xy_0': ra_at_xy_0, + kwargs_data = {'image_data': np.zeros((self.numpix, self.numpix)), + 'ra_at_xy_0': ra_at_xy_0, 'dec_at_xy_0': dec_at_xy_0, 'transform_pix2angle': transform_pix2angle, 'background_rms': self.background_noise, 'exposure_time': scaled_exposure_time} - return kwargs_data + return kwargs_data \ No newline at end of file diff --git a/lenstronomy/SimulationAPI/model_api.py b/lenstronomy/SimulationAPI/model_api.py index 49bb1e883..6f156d22f 100644 --- a/lenstronomy/SimulationAPI/model_api.py +++ b/lenstronomy/SimulationAPI/model_api.py @@ -23,20 +23,20 @@ def __init__(self, lens_model_list=None, z_lens=None, z_source=None, lens_redshi :param lens_model_list: list of strings with lens model names :param z_lens: redshift of the deflector (only considered when operating in single plane mode). - Is only needed for specific functions that require a cosmology. + Is only needed for specific functions that require a cosmology. :param z_source: redshift of the source: Needed in multi_plane option only, - not required for the core functionalities in the single plane mode. This will be the redshift of the source - plane (if not further specified the 'source_redshift_list') and the point source redshift - (regardless of 'source_redshift_list') + not required for the core functionalities in the single plane mode. This will be the redshift of the source + plane (if not further specified the 'source_redshift_list') and the point source redshift + (regardless of 'source_redshift_list') :param lens_redshift_list: list of deflector redshift (corresponding to the lens model list), - only applicable in multi_plane mode. + only applicable in multi_plane mode. :param source_light_model_list: list of strings with source light model names (lensed light profiles) :param lens_light_model_list: list of strings with lens light model names (not lensed light profiles) :param point_source_model_list: list of strings with point source model names :param source_redshift_list: list of redshifts of the source profiles (optional) :param cosmo: instance of the astropy cosmology class. If not specified, uses the default cosmology. :param z_source_convention: float, redshift of a source to define the reduced deflection angles of the lens - models. If None, 'z_source' is used. + models. If None, 'z_source' is used. """ if lens_model_list is None: lens_model_list = [] diff --git a/lenstronomy/SimulationAPI/observation_api.py b/lenstronomy/SimulationAPI/observation_api.py index bb2630baf..7f89c527f 100644 --- a/lenstronomy/SimulationAPI/observation_api.py +++ b/lenstronomy/SimulationAPI/observation_api.py @@ -43,9 +43,9 @@ def __init__(self, exposure_time, sky_brightness=None, seeing=None, num_exposure :param num_exposures: number of exposures that are combined :param psf_type: string, type of PSF ('GAUSSIAN' and 'PIXEL' supported) :param kernel_point_source: 2d numpy array, model of PSF centered with odd number of pixels per axis - (optional when psf_type='PIXEL' is chosen) + (optional when psf_type='PIXEL' is chosen) :param point_source_supersampling_factor: int, supersampling factor of kernel_point_source - (optional when psf_type='PIXEL' is chosen) + (optional when psf_type='PIXEL' is chosen) """ self._exposure_time = exposure_time self._sky_brightness_ = sky_brightness @@ -155,7 +155,7 @@ def __init__(self, pixel_scale, exposure_time, magnitude_zero_point, read_noise= :param magnitude_zero_point: magnitude in which 1 count (e-) per second per arcsecond square is registered :param num_exposures: number of exposures that are combined :param point_source_supersampling_factor: int, supersampling factor of kernel_point_source - (optional when psf_type='PIXEL' is chosen) + (optional when psf_type='PIXEL' is chosen) :param data_count_unit: string, unit of the data (not noise properties - see other definitions), 'e-': (electrons assumed to be IID), 'ADU': (analog-to-digital unit) diff --git a/lenstronomy/Util/image_util.py b/lenstronomy/Util/image_util.py index 65c160277..981e0cd3e 100644 --- a/lenstronomy/Util/image_util.py +++ b/lenstronomy/Util/image_util.py @@ -3,7 +3,6 @@ import numpy as np from scipy import ndimage from scipy import interpolate -from scipy.ndimage import interpolation as interp import copy import lenstronomy.Util.util as util @@ -15,6 +14,7 @@ def add_layer2image(grid2d, x_pos, y_pos, kernel, order=1): """ adds a kernel on the grid2d image at position x_pos, y_pos with an interpolated subgrid pixel shift of order=order + :param grid2d: 2d pixel grid (i.e. image) :param x_pos: x-position center (pixel coordinate) of the layer to be added :param y_pos: y-position center (pixel coordinate) of the layer to be added @@ -27,7 +27,7 @@ def add_layer2image(grid2d, x_pos, y_pos, kernel, order=1): y_int = int(round(y_pos)) shift_x = x_int - x_pos shift_y = y_int - y_pos - kernel_shifted = interp.shift(kernel, shift=[-shift_y, -shift_x], order=order) + kernel_shifted = ndimage.shift(kernel, shift=[-shift_y, -shift_x], order=order) return add_layer2image_int(grid2d, x_int, y_int, kernel_shifted) @@ -35,6 +35,7 @@ def add_layer2image(grid2d, x_pos, y_pos, kernel, order=1): def add_layer2image_int(grid2d, x_pos, y_pos, kernel): """ adds a kernel on the grid2d image at position x_pos, y_pos at integer positions of pixel + :param grid2d: 2d pixel grid (i.e. image) :param x_pos: x-position center (pixel coordinate) of the layer to be added :param y_pos: y-position center (pixel coordinate) of the layer to be added @@ -121,6 +122,7 @@ def rotateImage(img, angle): def re_size_array(x_in, y_in, input_values, x_out, y_out): """ resizes 2d array (i.e. image) to new coordinates. So far only works with square output aligned with coordinate axis. + :param x_in: :param y_in: :param input_values: @@ -138,6 +140,7 @@ def re_size_array(x_in, y_in, input_values, x_out, y_out): def symmetry_average(image, symmetry): """ symmetry averaged image + :param image: :param symmetry: :return: @@ -176,8 +179,6 @@ def coordInImage(x_coord, y_coord, num_pix, deltapix): checks whether image positions are within the pixel image in units of arcsec if not: remove it - :param imcoord: image coordinate (in units of angels) [[x,y,delta,magnification][...]] - :type imcoord: (n,4) numpy array :returns: image positions within the pixel image """ idex = [] @@ -243,14 +244,7 @@ def rebin_image(bin_size, image, wht_map, sigma_bkg, ra_coords, dec_coords, idex def rebin_coord_transform(factor, x_at_radec_0, y_at_radec_0, Mpix2coord, Mcoord2pix): """ adopt coordinate system and transformation between angular and pixel coordinates of a re-binned image - :param bin_size: - :param ra_0: - :param dec_0: - :param x_0: - :param y_0: - :param Matrix: - :param Matrix_inv: - :return: + """ factor = int(factor) Mcoord2pix_resized = Mcoord2pix / factor @@ -265,7 +259,7 @@ def rebin_coord_transform(factor, x_at_radec_0, y_at_radec_0, Mpix2coord, Mcoord def stack_images(image_list, wht_list, sigma_list): """ stacks images and saves new image as a fits file - :param image_name_list: list of image_names to be stacked + :return: """ image_stacked = np.zeros_like(image_list[0]) @@ -286,6 +280,7 @@ def cut_edges(image, num_pix): """ cuts out the edges of a 2d image and returns re-sized image to numPix center is well defined for odd pixel sizes. + :param image: 2d numpy array :param num_pix: square size of cut out image :return: cutout image with size numPix diff --git a/lenstronomy/Util/kernel_util.py b/lenstronomy/Util/kernel_util.py index 90f4f8ed2..fb8f91422 100644 --- a/lenstronomy/Util/kernel_util.py +++ b/lenstronomy/Util/kernel_util.py @@ -4,6 +4,7 @@ import numpy as np import copy import scipy.ndimage.interpolation as interp +from scipy import ndimage import lenstronomy.Util.util as util import lenstronomy.Util.image_util as image_util from lenstronomy.LightModel.Profiles.gaussian import Gaussian @@ -19,8 +20,8 @@ def de_shift_kernel(kernel, shift_x, shift_y, iterations=20, fractional_step_siz de-shifts a shifted kernel to the center of a pixel. This is performed iteratively. The input kernel is the solution of a linear interpolated shift of a sharper kernel centered in the middle of the - pixel. To find the de-shifted kernel, we perform an iterative correction of proposed de-shifted kernels and compare - its shifted version with the input kernel. + pixel. To find the de-shifted kernel, we perform an iterative correction of proposed de-shifted kernels and compare + its shifted version with the input kernel. :param kernel: (shifted) kernel, e.g. a star in an image that is not centered in the pixel grid :param shift_x: x-offset relative to the center of the pixel (sub-pixel shift) @@ -37,11 +38,11 @@ def de_shift_kernel(kernel, shift_x, shift_y, iterations=20, fractional_step_siz int_shift_y = int(round(shift_y)) frac_y_shift = shift_y - int_shift_y kernel_init = copy.deepcopy(kernel_new) - kernel_init_shifted = copy.deepcopy(interp.shift(kernel_init, shift=[int_shift_y, int_shift_x], order=1)) - kernel_new = interp.shift(kernel_new, shift=[int_shift_y, int_shift_x], order=1) + kernel_init_shifted = copy.deepcopy(ndimage.shift(kernel_init, shift=[int_shift_y, int_shift_x], order=1)) + kernel_new = ndimage.shift(kernel_new, shift=[int_shift_y, int_shift_x], order=1) norm = np.sum(kernel_init_shifted) for i in range(iterations): - kernel_shifted_inv = interp.shift(kernel_new, shift=[-frac_y_shift, -frac_x_shift], order=1) + kernel_shifted_inv = ndimage.shift(kernel_new, shift=[-frac_y_shift, -frac_x_shift], order=1) delta = kernel_init_shifted - kernel_norm(kernel_shifted_inv) * norm kernel_new += delta * fractional_step_size kernel_new = kernel_norm(kernel_new) * norm @@ -149,6 +150,7 @@ def subgrid_kernel(kernel, subgrid_res, odd=False, num_iter=100): def kernel_pixelsize_change(kernel, deltaPix_in, deltaPix_out): """ change the pixel size of a given kernel + :param kernel: :param deltaPix_in: :param deltaPix_out: @@ -169,6 +171,7 @@ def kernel_pixelsize_change(kernel, deltaPix_in, deltaPix_out): def cut_psf(psf_data, psf_size): """ cut the psf properly + :param psf_data: image of PSF :param psf_size: size of psf :return: re-sized and re-normalized PSF @@ -371,6 +374,7 @@ def averaging_even_kernel(kernel_high_res, subgrid_res): def cutout_source(x_pos, y_pos, image, kernelsize, shift=True): """ cuts out point source (e.g. PSF estimate) out of image and shift it to the center of a pixel + :param x_pos: :param y_pos: :param image: @@ -434,6 +438,7 @@ def fwhm_kernel(kernel): def estimate_amp(data, x_pos, y_pos, psf_kernel): """ estimates the amplitude of a point source located at x_pos, y_pos + :param data: :param x_pos: :param y_pos: diff --git a/lenstronomy/Util/prob_density.py b/lenstronomy/Util/prob_density.py index 518a43ec2..ad662c46f 100644 --- a/lenstronomy/Util/prob_density.py +++ b/lenstronomy/Util/prob_density.py @@ -16,6 +16,7 @@ def pdf(self, x, e=0., w=1., a=0.): """ probability density function see: https://en.wikipedia.org/wiki/Skew_normal_distribution + :param x: input value :param e: :param w: @@ -28,6 +29,7 @@ def pdf(self, x, e=0., w=1., a=0.): def pdf_skew(self, x, mu, sigma, skw): """ function with different parameterisation + :param x: :param mu: mean :param sigma: sigma @@ -61,6 +63,7 @@ def _alpha_delta(self, delta): def _w_sigma_delta(self, sigma, delta): """ invert variance + :param sigma: :param delta: :return: w parameter @@ -84,6 +87,7 @@ def _e_mu_w_delta(self, mu, w, delta): def map_mu_sigma_skw(self, mu, sigma, skw): """ map to parameters e, w, a + :param mu: mean :param sigma: standard deviation :param skw: skewness @@ -125,6 +129,7 @@ def compute_lower_upper_errors(sample, num_sigma=1): """ computes the upper and lower sigma from the median value. This functions gives good error estimates for skewed pdf's + :param sample: 1-D sample :param num_sigma: integer, number of sigmas to be returned :return: median, lower_sigma, upper_sigma diff --git a/lenstronomy/Util/util.py b/lenstronomy/Util/util.py index 009010e18..e1d7b1fbb 100644 --- a/lenstronomy/Util/util.py +++ b/lenstronomy/Util/util.py @@ -95,6 +95,7 @@ def map_coord2pix(ra, dec, x_0, y_0, M): """ this routines performs a linear transformation between two coordinate systems. Mainly used to transform angular into pixel coordinates in an image + :param ra: ra coordinates :param dec: dec coordinates :param x_0: pixel value in x-axis of ra,dec = 0,0 @@ -226,6 +227,7 @@ def make_grid(numPix, deltapix, subgrid_res=1, left_lower=False): def make_grid_transformed(numPix, Mpix2Angle): """ returns grid with linear transformation (deltaPix and rotation) + :param numPix: number of Pixels :param Mpix2Angle: 2-by-2 matrix to mat a pixel to a coordinate :return: coordinate grid @@ -281,7 +283,6 @@ def grid_from_coordinate_transform(nx, ny, Mpix2coord, ra_at_xy_0, dec_at_xy_0): """ return a grid in x and y coordinates that satisfy the coordinate system - :param nx: number of pixels in x-axis :param ny: number of pixels in y-axis :param Mpix2coord: transformation matrix (2x2) of pixels into coordinate displacements @@ -303,6 +304,7 @@ def grid_from_coordinate_transform(nx, ny, Mpix2coord, ra_at_xy_0, dec_at_xy_0): def get_axes(x, y): """ computes the axis x and y of a given 2d grid + :param x: :param y: :return: @@ -321,6 +323,7 @@ def get_axes(x, y): def averaging(grid, numGrid, numPix): """ resize 2d pixel grid with numGrid to numPix and averages over the pixels + :param grid: higher resolution pixel grid :param numGrid: number of pixels per axis in the high resolution input image :param numPix: lower number of pixels per axis in the output image (numGrid/numPix is integer number) @@ -403,6 +406,7 @@ def compare_distance(x_mapped, y_mapped): def min_square_dist(x_1, y_1, x_2, y_2): """ return minimum of quadratic distance of pairs (x1, y1) to pairs (x2, y2) + :param x_1: :param y_1: :param x_2: @@ -467,6 +471,7 @@ def select_best(array, criteria, num_select, highest=True): def points_on_circle(radius, num_points, connect_ends=True): """ returns a set of uniform points around a circle + :param radius: radius of the circle :param num_points: number of points on the circle :param connect_ends: boolean, if True, start and end point are the same @@ -480,6 +485,40 @@ def points_on_circle(radius, num_points, connect_ends=True): y_coord = np.sin(angle) * radius return x_coord, y_coord +@export +@jit() +def local_minima_2d(a, x, y): + """ + finds (local) minima in a 2d grid + applies less rigid criteria for maximum without second-order tangential minima criteria + + :param a: 1d array of displacements from the source positions + :type a: numpy array with length numPix**2 in float + :param x: 1d coordinate grid in x-direction + :type x: numpy array with length numPix**2 in float + :param y: 1d coordinate grid in x-direction + :type y: numpy array with length numPix**2 in float + :returns: array of indices of local minima, values of those minima + :raises: AttributeError, KeyError + """ + dim = int(np.sqrt(len(a))) + values = [] + x_mins = [] + y_mins = [] + for i in range(dim + 1, len(a) - dim - 1): + if (a[i] < a[i - 1] + and a[i] < a[i + 1] + and a[i] < a[i - dim] + and a[i] < a[i + dim] + and a[i] < a[i - (dim - 1)] + and a[i] < a[i - (dim + 1)] + and a[i] < a[i + (dim - 1)] + and a[i] < a[i + (dim + 1)]): + x_mins.append(x[i]) + y_mins.append(y[i]) + values.append(a[i]) + return np.array(x_mins), np.array(y_mins), np.array(values) + @export @jit() @@ -540,7 +579,7 @@ def neighborSelect(a, x, y): def fwhm2sigma(fwhm): """ - :param fwhm: full-widt-half-max value + :param fwhm: full-width-half-max value :return: gaussian sigma (sqrt(var)) """ sigma = fwhm / (2 * np.sqrt(2 * np.log(2))) @@ -585,6 +624,7 @@ def hyper2F2_array(a, b, c, d, x): def make_subgrid(ra_coord, dec_coord, subgrid_res=2): """ return a grid with subgrid resolution + :param ra_coord: :param dec_coord: :param subgrid_res: @@ -620,6 +660,7 @@ def make_subgrid(ra_coord, dec_coord, subgrid_res=2): def convert_bool_list(n, k=None): """ returns a bool list of the length of the lens models + if k = None: returns bool list with True's if k is int, returns bool list with False's but k'th is True if k is a list of int, e.g. [0, 3, 5], returns a bool list with True's in the integers listed and False elsewhere diff --git a/lenstronomy/Workflow/fitting_sequence.py b/lenstronomy/Workflow/fitting_sequence.py index 0908232e4..a51ef6652 100644 --- a/lenstronomy/Workflow/fitting_sequence.py +++ b/lenstronomy/Workflow/fitting_sequence.py @@ -40,7 +40,7 @@ def __init__(self, kwargs_data_joint, kwargs_model, kwargs_constraints, kwargs_l 'extinction_model': [kwargs_init, kwargs_sigma, kwargs_fixed, kwargs_lower, kwargs_upper] 'special': [kwargs_init, kwargs_sigma, kwargs_fixed, kwargs_lower, kwargs_upper] :param mpi: MPI option (bool), if True, will launch an MPI Pool job for the steps in the fitting sequence where - possible + possible :param verbose: bool, if True prints temporary results and indicators of the fitting process """ self.kwargs_data_joint = kwargs_data_joint @@ -164,7 +164,8 @@ def best_fit_likelihood(self): @property def bic(self): """ - returns the bayesian information criterion of the model. + Bayesian information criterion (BIC) of the model. + :return: bic value, float """ num_data = self.likelihoodModule.num_data @@ -493,6 +494,7 @@ def update_settings(self, kwargs_model=None, kwargs_constraints=None, kwargs_lik def set_param_value(self, **kwargs): """ Set a parameter to a specific value. `kwargs` are below. + :param lens: [[i_model, ['param1', 'param2',...], [...]] :type lens: :param source: [[i_model, ['param1', 'param2',...], [...]] diff --git a/lenstronomy/Workflow/psf_fitting.py b/lenstronomy/Workflow/psf_fitting.py index 343a3f6f0..af55b8b85 100644 --- a/lenstronomy/Workflow/psf_fitting.py +++ b/lenstronomy/Workflow/psf_fitting.py @@ -7,6 +7,7 @@ import numpy as np import copy import scipy.ndimage.interpolation as interp +from scipy import ndimage __all__ = ['PsfFitting'] @@ -315,11 +316,11 @@ def psf_estimate_individual(self, ra_image, dec_image, point_amp, residuals, cut shift_y = (y_int - y_[l]) * supersampling_factor # for odd number super-sampling if supersampling_factor % 2 == 1: - residuals_shifted = interp.shift(residual_cutout_mask, shift=[shift_y, shift_x], order=1) + residuals_shifted = ndimage.shift(residual_cutout_mask, shift=[shift_y, shift_x], order=1) else: # for even number super-sampling half a super-sampled pixel offset needs to be performed - residuals_shifted = interp.shift(residual_cutout_mask, shift=[shift_y - 0.5, shift_x - 0.5], order=1) + residuals_shifted = ndimage.shift(residual_cutout_mask, shift=[shift_y - 0.5, shift_x - 0.5], order=1) # and the last column and row need to be removed residuals_shifted = residuals_shifted[:-1, :-1] @@ -379,7 +380,7 @@ def cutout_psf_single(x, y, image, mask, kernel_size, kernel_init): # shift the initial kernel to the shift of the star shift_x = x_int - x shift_y = y_int - y - kernel_shifted = interp.shift(kernel_enlarged, shift=[-shift_y, -shift_x], order=1) + kernel_shifted = ndimage.shift(kernel_enlarged, shift=[-shift_y, -shift_x], order=1) # compute normalization of masked and unmasked region of the shifted kernel # norm_masked = np.sum(kernel_shifted[mask_i == 0]) norm_unmasked = np.sum(kernel_shifted[mask_cutout == 1]) @@ -401,10 +402,11 @@ def cutout_psf_single(x, y, image, mask, kernel_size, kernel_init): def combine_psf(kernel_list_new, kernel_old, factor=1., stacking_option='median', symmetry=1): """ updates psf estimate based on old kernel and several new estimates + :param kernel_list_new: list of new PSF kernels estimated from the point sources in the image (un-normalized) :param kernel_old: old PSF kernel :param factor: weight of updated estimate based on new and old estimate, factor=1 means new estimate, - factor=0 means old estimate + factor=0 means old estimate :param stacking_option: option of stacking, mean or median :param symmetry: imposed symmetry of PSF estimate :return: updated PSF estimate and error_map associated with it @@ -514,7 +516,7 @@ def error_map_estimate(self, kernel, star_cutout_list, amp, x_pos, y_pos, error_ y_int = int(round(y)) shift_x = x_int - x shift_y = y_int - y - kernel_shifted = interp.shift(kernel, shift=[-shift_y, -shift_x], order=1) + kernel_shifted = ndimage.shift(kernel, shift=[-shift_y, -shift_x], order=1) # multiply kernel with amplitude model = kernel_shifted * amp_i # compute residuals diff --git a/lenstronomy/Workflow/update_manager.py b/lenstronomy/Workflow/update_manager.py index 2234bf0e8..1dd3a8771 100644 --- a/lenstronomy/Workflow/update_manager.py +++ b/lenstronomy/Workflow/update_manager.py @@ -300,9 +300,7 @@ def update_fixed(self, lens_add_fixed=None, source_add_fixed=None, lens_light_ad if special_add_fixed is None: special_add_fixed = [] for param_name in special_add_fixed: - if param_name in special_fixed: - pass - else: + if param_name not in special_fixed: special_fixed[param_name] = special_temp[param_name] if special_remove_fixed is None: special_remove_fixed = [] diff --git a/requirements.txt b/requirements.txt index 883bb3ed3..4b67d92c7 100644 --- a/requirements.txt +++ b/requirements.txt @@ -15,7 +15,7 @@ pyyaml pyxdg h5py zeus-mcmc -nautilus-sampler==0.1.0 +nautilus-sampler>=0.2.1 schwimmbad diff --git a/setup.py b/setup.py index e547a199c..79ff446c3 100644 --- a/setup.py +++ b/setup.py @@ -56,7 +56,7 @@ def run_tests(self): ] tests_require = ['pytest>=2.3', "mock", 'colossus==1.3.0', 'slitronomy==0.3.2', 'emcee>=3.0.0', 'dynesty', 'nestcheck', 'pymultinest', 'zeus-mcmc>=2.4.0', - 'nautilus-sampler==0.1.0', + 'nautilus-sampler>=0.2.1', ] PACKAGE_PATH = os.path.abspath(os.path.join(__file__, os.pardir)) diff --git a/test/test_Cosmo/test_cosmo_solver.py b/test/test_Cosmo/test_cosmo_solver.py index e9d54a37d..1d3c313e1 100644 --- a/test/test_Cosmo/test_cosmo_solver.py +++ b/test/test_Cosmo/test_cosmo_solver.py @@ -74,6 +74,7 @@ def setup(self): self.z_d, self.z_s = 0.295, 0.658 self.invertCosmo = InvertCosmo(z_d=self.z_d, z_s=self.z_s, H0_range=np.linspace(10, 100, 50), omega_m_range=np.linspace(0.05, 1, 50)) + self.invertCosmo_default = InvertCosmo(z_d=self.z_d, z_s=self.z_s) def test_get_cosmo(self): H0 = 80 @@ -83,6 +84,10 @@ def test_get_cosmo(self): npt.assert_almost_equal(H0_new, H0, decimal=1) npt.assert_almost_equal(omega_m_new, omega_m, decimal=3) + H0_new, omega_m_new = self.invertCosmo_default.get_cosmo(Dd, Ds_Dds) + npt.assert_almost_equal(H0_new, H0, decimal=1) + npt.assert_almost_equal(omega_m_new, omega_m, decimal=3) + H0_new, omega_m_new = self.invertCosmo.get_cosmo(Dd=1, Ds_Dds=1) assert H0_new == -1 diff --git a/test/test_Cosmo/test_lens_cosmo.py b/test/test_Cosmo/test_lens_cosmo.py index 83b5d1143..db4fc6634 100644 --- a/test/test_Cosmo/test_lens_cosmo.py +++ b/test/test_Cosmo/test_lens_cosmo.py @@ -108,6 +108,16 @@ def test_a_z(self): a = self.lensCosmo.a_z(z=1) npt.assert_almost_equal(a, 0.5) + def test_sersic_m_star2k_eff(self): + m_star = 10**11.5 + R_sersic = 1 + n_sersic = 4 + k_eff = self.lensCosmo.sersic_m_star2k_eff(m_star, R_sersic, n_sersic) + npt.assert_almost_equal(k_eff, 0.1294327891669961, decimal=5) + + m_star_out = self.lensCosmo.sersic_k_eff2m_star(k_eff, R_sersic, n_sersic) + npt.assert_almost_equal(m_star_out, m_star, decimal=6) + if __name__ == '__main__': pytest.main() diff --git a/test/test_Data/test_psf.py b/test/test_Data/test_psf.py index 40e19dd58..1084f4779 100644 --- a/test/test_Data/test_psf.py +++ b/test/test_Data/test_psf.py @@ -31,8 +31,8 @@ def test_kernel_point_source(self): kernel_point_source = psf_class.kernel_point_source assert len(kernel_point_source) == 13 kernel_super = psf_class.kernel_point_source_supersampled(supersampling_factor=3) - assert np.sum(kernel_point_source) == np.sum(kernel_super) - assert np.sum(kernel_point_source) == 1 + npt.assert_almost_equal(np.sum(kernel_point_source), np.sum(kernel_super), decimal=9) + npt.assert_almost_equal(np.sum(kernel_point_source), 1, decimal=9) def test_kernel_subsampled(self): deltaPix = 0.05 # pixel size of image diff --git a/test/test_LensModel/test_Profiles/test_cored_steep_ellipsoid.py b/test/test_LensModel/test_Profiles/test_cored_steep_ellipsoid.py index 13d65d164..fdbecf81e 100644 --- a/test/test_LensModel/test_Profiles/test_cored_steep_ellipsoid.py +++ b/test/test_LensModel/test_Profiles/test_cored_steep_ellipsoid.py @@ -4,6 +4,7 @@ import numpy as np import pytest import numpy.testing as npt +from lenstronomy.Util import param_util class TestCSP(object): @@ -12,7 +13,7 @@ class TestCSP(object): """ def setup(self): from lenstronomy.LensModel.Profiles.cored_steep_ellipsoid import CSE - self.CSP = CSE() + self.CSP = CSE(axis='product_avg') def test_function(self): @@ -43,6 +44,33 @@ def test_hessian(self): npt.assert_almost_equal(f_xy, [-0.13493, -0.], decimal=5) npt.assert_almost_equal(f_yy, [-0.03315, 0.27639], decimal=5) + def test_ellipticity(self): + """ + test the definition of the ellipticity normalization (along major axis or product averaged axes) + """ + x, y = np.linspace(start=0.001, stop=10, num=100), np.zeros(100) + kwargs_round = {'a': 2, 's': 1, 'e1': 0., 'e2': 0., 'center_x': 0, 'center_y': 0} + phi_q, q = param_util.ellipticity2phi_q(0.3, 0) + kwargs = {'a': 2, 's': 1, 'e1': 0.3, 'e2': 0., 'center_x': 0, 'center_y': 0} + + f_xx, f_xy, f_yx, f_yy = self.CSP.hessian(x, y, **kwargs_round) + kappa_round = 1. / 2 * (f_xx + f_yy) + + f_xx, f_xy, f_yx, f_yy = self.CSP.hessian(x, y, **kwargs) + kappa_major = 1. / 2 * (f_xx + f_yy) + + f_xx, f_xy, f_yx, f_yy = self.CSP.hessian(y, x, **kwargs) + kappa_minor = 1. / 2 * (f_xx + f_yy) + + # import matplotlib.pyplot as plt + # plt.plot(x, kappa_major/kappa_round, ',-', label='major/round', alpha=0.5) + # plt.plot(x, kappa_minor/kappa_round, '--', label='minor/round', alpha=0.5) + # + # plt.plot(x, np.sqrt(kappa_minor*kappa_major)/kappa_round,label='prod/kappa_round') + # plt.legend() + # plt.show() + + npt.assert_almost_equal(kappa_round,np.sqrt(kappa_minor*kappa_major), decimal=1) if __name__ == '__main__': pytest.main() diff --git a/test/test_LensModel/test_Profiles/test_gaussian.py b/test/test_LensModel/test_Profiles/test_gaussian.py index 1464dcb4c..0efcb56d6 100644 --- a/test/test_LensModel/test_Profiles/test_gaussian.py +++ b/test/test_LensModel/test_Profiles/test_gaussian.py @@ -68,9 +68,9 @@ def test_hessian(self): assert values[0][0] == 0. assert values[3][0] == -np.exp(-1./2) assert values[1][0] == 0. - assert values[0][1] == 0.40600584970983811 - assert values[3][1] == -0.1353352832366127 - assert values[1][1] == 0. + npt.assert_almost_equal(values[0][1], 0.40600584970983811, decimal=9) + npt.assert_almost_equal(values[3][1], -0.1353352832366127, decimal=9) + npt.assert_almost_equal(values[1][1], 0, decimal=9) class TestGaussianKappa(object): diff --git a/test/test_LensModel/test_Profiles/test_nfw_ellipse_cse.py b/test/test_LensModel/test_Profiles/test_nfw_ellipse_cse.py index fd107f6e5..f8749f888 100644 --- a/test/test_LensModel/test_Profiles/test_nfw_ellipse_cse.py +++ b/test/test_LensModel/test_Profiles/test_nfw_ellipse_cse.py @@ -2,6 +2,8 @@ from lenstronomy.LensModel.Profiles.nfw import NFW from lenstronomy.LensModel.Profiles.nfw_ellipse_cse import NFW_ELLIPSE_CSE +from lenstronomy.Analysis.lens_profile import LensProfileAnalysis +from lenstronomy.LensModel.lens_model import LensModel import numpy as np import numpy.testing as npt @@ -67,6 +69,45 @@ def test_mass_3d_lens(self): m_3d_cse = self.nfw_cse.mass_3d_lens(R, Rs, alpha_Rs) npt.assert_almost_equal(m_3d_nfw, m_3d_cse, decimal=8) + def test_ellipticity(self): + """ + test the definition of the ellipticity normalization (along major axis or product averaged axes) + """ + x, y = np.linspace(start=0.001, stop=10, num=100), np.zeros(100) + kwargs_round = {'alpha_Rs': 0.5, 'Rs': 2, 'center_x': 0, 'center_y': 0, 'e1': 0, 'e2': 0} + kwargs = {'alpha_Rs': 0.5, 'Rs': 2, 'center_x': 0, 'center_y': 0, 'e1': 0.3, 'e2': 0} + + f_xx, f_xy, f_yx, f_yy = self.nfw_cse.hessian(x, y, **kwargs_round) + kappa_round = 1. / 2 * (f_xx + f_yy) + + f_xx, f_xy, f_yx, f_yy = self.nfw_cse.hessian(x, y, **kwargs) + kappa_major = 1. / 2 * (f_xx + f_yy) + + f_xx, f_xy, f_yx, f_yy = self.nfw_cse.hessian(y, x, **kwargs) + kappa_minor = 1. / 2 * (f_xx + f_yy) + + npt.assert_almost_equal(np.sqrt(kappa_minor * kappa_major),kappa_round, decimal=2) + + # import matplotlib.pyplot as plt + # plt.plot(x, kappa_round/kappa_round, ':', label='round', alpha=0.5) + # plt.plot(x, kappa_major/kappa_round, ',-', label='major', alpha=0.5) + # plt.plot(x, kappa_minor/kappa_round, '--', label='minor', alpha=0.5) + # plt.plot(x, np.sqrt(kappa_minor * kappa_major)/kappa_round, '--', label='prod', alpha=0.5) + # plt.plot(x, np.sqrt(kappa_minor**2 + kappa_major**2) / kappa_round / 2, '--', label='square', alpha=0.5) + # plt.legend() + # plt.show() + def test_einstein_rad(self): + """ + test that the Einstein radius doesn't change significantly with ellipticity + """ + kwargs_round = {'alpha_Rs': 0.5, 'Rs': 2, 'center_x': 0, 'center_y': 0, 'e1': 0, 'e2': 0} + kwargs = {'alpha_Rs': 0.5, 'Rs': 2, 'center_x': 0, 'center_y': 0, 'e1': 0.3, 'e2': 0} + LensMod=LensModel(['NFW_ELLIPSE_CSE']) + LensAn=LensProfileAnalysis(LensMod) + r_Ein_round=LensAn.effective_einstein_radius([kwargs_round]) + r_Ein_ell=LensAn.effective_einstein_radius([kwargs]) + npt.assert_almost_equal(r_Ein_round,r_Ein_ell,decimal=1) + if __name__ == '__main__': pytest.main() diff --git a/test/test_LensModel/test_Profiles/test_shapelet_pot_cartesian.py b/test/test_LensModel/test_Profiles/test_shapelet_pot_cartesian.py index cdb8f16ff..a75f109c6 100644 --- a/test/test_LensModel/test_Profiles/test_shapelet_pot_cartesian.py +++ b/test/test_LensModel/test_Profiles/test_shapelet_pot_cartesian.py @@ -22,14 +22,14 @@ def test_function(self): beta = 1. coeffs = (1., 1.) values = self.cartShapelets.function(x, y, coeffs, beta) - assert values[0] == 0.11180585426466891 + npt.assert_almost_equal(values[0], 0.11180585426466888, decimal=9) x = 1. y = 2. beta = 1. coeffs = (1., 1.) values = self.cartShapelets.function(x, y, coeffs, beta) - assert values == 0.11180585426466891 + npt.assert_almost_equal(values, 0.11180585426466891, decimal=9) x = np.array([0]) y = np.array([0]) @@ -40,11 +40,11 @@ def test_function(self): coeffs = (1, 1., 0, 0, 1, 1) values = self.cartShapelets.function(x, y, coeffs, beta) - assert values[0] == 0.16524730314632363 + npt.assert_almost_equal(values[0], 0.16524730314632363, decimal=9) coeffs = (1, 1., 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0) values = self.cartShapelets.function(x, y, coeffs, beta) - assert values[0] == 0.16524730314632363 + npt.assert_almost_equal(values[0], 0.16524730314632363, decimal=9) coeffs = (0., 0., 0, 0, 0., 0., 0, 0, 0, 0, 0, 0, 0, 0, 0) values = self.cartShapelets.function(x, y, coeffs, beta) diff --git a/test/test_LensModel/test_Solver/test_lens_equation_solver.py b/test/test_LensModel/test_Solver/test_lens_equation_solver.py index 90c677f50..b45b88e68 100644 --- a/test/test_LensModel/test_Solver/test_lens_equation_solver.py +++ b/test/test_LensModel/test_Solver/test_lens_equation_solver.py @@ -145,6 +145,21 @@ def test_analytical_lens_equation_solver(self): for x, y in zip(x_pos_ls, y_pos_ls): # Check if it found all solutions lenstronomy found assert np.sqrt((x-x_pos)**2+(y-y_pos)**2).min() < 1e-8 + # here we test with shear and mass profile centroids not aligned + lensModel = LensModel(['EPL_NUMBA', 'SHEAR']) + lensEquationSolver = LensEquationSolver(lensModel) + sourcePos_x = 0.03 + sourcePos_y = 0.0 + kwargs_lens = [{'theta_E': 1., 'gamma': 2.2, 'center_x': 0.01, 'center_y': 0.02, 'e1': 0.01, 'e2': 0.05}, + {'gamma1': -0.04, 'gamma2': -0.1, 'ra_0': 0.0, 'dec_0': 1.0}] + + x_pos, y_pos = lensEquationSolver.image_position_from_source(sourcePos_x, sourcePos_y, kwargs_lens, + solver='analytical') + source_x, source_y = lensModel.ray_shooting(x_pos, y_pos, kwargs_lens) + assert len(source_x) == len(source_y) >= 2 + npt.assert_almost_equal(sourcePos_x, source_x, decimal=10) + npt.assert_almost_equal(sourcePos_y, source_y, decimal=10) + def test_caustics(self): lm = LensModel(['EPL_NUMBA', 'SHEAR']) leqs = LensEquationSolver(lm) @@ -159,7 +174,8 @@ def test_caustics(self): lensplane_cut = caustics_epl_shear(kwargs, return_which='cut', sourceplane=False) twoimg = caustics_epl_shear(kwargs, return_which='double') fourimg = caustics_epl_shear(kwargs, return_which='quad') - assert np.abs(lm.magnification(*lensplane_caus, kwargs)).min() > 1e12 + min_mag = np.abs(lm.magnification(*lensplane_caus, kwargs)).min() + assert min_mag > 1e12 assert np.abs(lm.magnification(*lensplane_cut, kwargs)).min() > 1e12 # Test whether the caustics indeed the number of images they say @@ -213,11 +229,6 @@ def test_assertions(self): with pytest.raises(ValueError): lensEquationSolver.image_position_from_source(0.1, 0., kwargs_lens, solver='nonexisting') - with pytest.raises(ValueError): - kwargs_lens[1]['ra_0']=0.1 - lensEquationSolver.image_position_from_source(0.1, 0., kwargs_lens, solver='analytical') - - if __name__ == '__main__': pytest.main() diff --git a/test/test_LensModel/test_single_plane.py b/test/test_LensModel/test_single_plane.py index e873eb492..915afe654 100644 --- a/test/test_LensModel/test_single_plane.py +++ b/test/test_LensModel/test_single_plane.py @@ -40,7 +40,7 @@ def test_mass_2d(self): lensModel = SinglePlane(['GAUSSIAN_KAPPA']) kwargs = [{'amp': 1., 'sigma': 2., 'center_x': 0., 'center_y': 0.}] output = lensModel.mass_2d(r=1, kwargs=kwargs) - assert output == 0.11750309741540453 + npt.assert_almost_equal(output, 0.11750309741540453, decimal=9) def test_density(self): theta_E = 1 diff --git a/test/test_LightModel/test_Profiles/test_interpolation.py b/test/test_LightModel/test_Profiles/test_interpolation.py index 265a99933..25f316491 100644 --- a/test/test_LightModel/test_Profiles/test_interpolation.py +++ b/test/test_LightModel/test_Profiles/test_interpolation.py @@ -3,6 +3,7 @@ from lenstronomy.LightModel.Profiles.interpolation import Interpol from lenstronomy.LightModel.Profiles.gaussian import Gaussian import lenstronomy.Util.util as util +import numpy as np class TestInterpol(object): @@ -55,6 +56,30 @@ def test_function(self): out_scaled = interp.function(x=2., y=0, image=image, scale=2, phi_G=0, center_x=0, center_y=0) assert out_scaled == out + def test_flux_normalization(self): + interp = Interpol() + delta_pix = 0.1 + len_x, len_y = 21, 21 + x, y = util.make_grid(numPix=(len_x, len_y), deltapix=delta_pix) + gauss = Gaussian() + flux = gauss.function(x, y, amp=1., center_x=0., center_y=0., sigma=0.3) + image = util.array2image(flux, nx=len_y, ny=len_x) + flux_total = np.sum(image) + + kwargs_interp = {'image': image, 'scale': delta_pix, 'phi_G': 0., 'center_x': 0., 'center_y': 0.} + image_interp = interp.function(x, y, **kwargs_interp) + flux_interp = np.sum(image_interp) + npt.assert_almost_equal(flux_interp, flux_total, decimal=3) + + # test for scale !=1 + # demands same surface brightness values. We rescale the pixel grid by the same amount as the image + scale = 0.5 + x, y = util.make_grid(numPix=(len_x, len_y), deltapix=delta_pix * scale) + kwargs_interp = {'image': image, 'scale': delta_pix * scale, 'phi_G': 0., 'center_x': 0., 'center_y': 0.} + output = interp.function(x, y, **kwargs_interp) + + npt.assert_almost_equal(output / image_interp, 1, decimal=5) + def test_delete_cache(self): x, y = util.make_grid(numPix=20, deltapix=1.) gauss = Gaussian() diff --git a/test/test_LightModel/test_Profiles/test_shapelets.py b/test/test_LightModel/test_Profiles/test_shapelets.py index eb07446a9..8e7e5aa9d 100644 --- a/test/test_LightModel/test_Profiles/test_shapelets.py +++ b/test/test_LightModel/test_Profiles/test_shapelets.py @@ -50,7 +50,7 @@ def test_shapelet_basis(self): beta = 1 numPix = 10 kernel_list = self.shapeletSet.shapelet_basis_2d(num_order, beta, numPix) - assert kernel_list[0][4, 4] == 0.43939128946772255 + npt.assert_almost_equal(kernel_list[0][4, 4], 0.4393912894677224, decimal=9) def test_decomposition(self): """ diff --git a/test/test_Sampling/test_Samplers/test_nautilus.py b/test/test_Sampling/test_Samplers/test_nautilus.py index 2694cb7a7..c15afce16 100644 --- a/test/test_Sampling/test_Samplers/test_nautilus.py +++ b/test/test_Sampling/test_Samplers/test_nautilus.py @@ -11,7 +11,7 @@ def import_fixture(simple_einstein_ring_likelihood_2d): """ - :param simple_einstein_ring_likelihood: fixture + :param simple_einstein_ring_likelihood_2d: fixture :return: """ likelihood, kwargs_truths = simple_einstein_ring_likelihood_2d @@ -46,6 +46,16 @@ def test_sampler(self, import_fixture): kwargs_run_fail['prior_type'] = 'wrong' assert_raises(ValueError, sampler.nautilus_sampling, **kwargs_run_fail) + def test_prior(self): + + num_param = 10 + from nautilus import Prior + prior = Prior() + + for i in range(num_param): + prior.add_parameter(dist=(0, 1)) + assert num_param == prior.dimensionality() + if __name__ == '__main__': pytest.main() diff --git a/test/test_Sampling/test_param_groups.py b/test/test_Sampling/test_param_groups.py new file mode 100644 index 000000000..3f164196d --- /dev/null +++ b/test/test_Sampling/test_param_groups.py @@ -0,0 +1,88 @@ +__author__ = 'jhodonnell' + + +import numpy as np +import numpy.testing as npt +import unittest +import pytest + +from lenstronomy.Sampling.param_group import ( + ModelParamGroup, SingleParam, ArrayParam + ) + + +class ExampleSingleParam(SingleParam): + param_names = ['sp1', 'sp2'] + _kwargs_lower = {'sp1': 0, 'sp2': 0} + _kwargs_upper = {'sp1': 10, 'sp2': 10} + + +class ExampleArrayParam(ArrayParam): + param_names = {'ap1': 1, 'ap2': 3} + _kwargs_lower = {'ap1': [0], 'ap2': [0]*3} + _kwargs_upper = {'ap1': [10], 'ap2': [10]*3} + + +class TestParamGroup(object): + def setup(self): + pass + + def test_single_param(self): + sp = ExampleSingleParam(on=True) + + num, names = sp.num_params({}) + assert num == 2 + assert names == ['sp1', 'sp2'] + + result = sp.set_params({'sp1': 2}, {'sp2': 3}) + assert result == [2] + + result = sp.set_params({'sp1': 2, 'sp2': 3}, {}) + assert result == [2, 3] + + kwargs, i = sp.get_params(result, i=0, kwargs_fixed={}) + assert kwargs['sp1'] == 2 + assert kwargs['sp2'] == 3 + + def test_array_param(self): + ap = ExampleArrayParam(on=True) + + num, names = ap.num_params({}) + assert num == 4 + assert names == ['ap1'] + ['ap2']*3 + + result = ap.set_params({'ap1': [2]}, {'ap2': [1, 1, 1]}) + assert result == [2] + + result = ap.set_params({'ap1': [2], 'ap2': [1, 2, 3]}, {}) + assert result == [2, 1, 2, 3] + + kwargs, i = ap.get_params(result, i=0, kwargs_fixed={}) + assert kwargs['ap1'] == [2] + assert kwargs['ap2'] == [1, 2, 3] + + def test_compose(self): + sp = ExampleSingleParam(on=True) + ap = ExampleArrayParam(on=True) + + num, names = ModelParamGroup.compose_num_params([sp, ap], kwargs_fixed={}) + assert num == 6 + assert names == sp.num_params({})[1] + ap.num_params({})[1] + + result = ModelParamGroup.compose_set_params( + [sp, ap], + { + 'sp1': 1, 'sp2': 2, + 'ap1': [3], 'ap2': [4, 5, 6] + }, + kwargs_fixed={} + ) + assert result == [1, 2, 3, 4, 5, 6] + + kwargs, i = ModelParamGroup.compose_get_params([sp, ap], result, i=0, kwargs_fixed={}) + assert kwargs['sp1'] == 1 + assert kwargs['ap2'] == [4, 5, 6] + + +if __name__ == '__main__': + pytest.main() diff --git a/test/test_Sampling/test_parameters.py b/test/test_Sampling/test_parameters.py index 107136856..0066a42d8 100644 --- a/test/test_Sampling/test_parameters.py +++ b/test/test_Sampling/test_parameters.py @@ -103,7 +103,7 @@ def test_get_cosmo(self): args = param_class.kwargs2args(kwargs_true_lens, kwargs_true_source, kwargs_lens_light=kwargs_true_lens_light, kwargs_ps=kwargs_true_ps, kwargs_special={'D_dt': 1000}) - assert param_class.specialParams._D_dt_sampling is True + assert param_class.specialParams._D_dt_sampling.on def test_mass_scaling(self): kwargs_model = {'lens_model_list': ['SIS', 'NFW', 'NFW', 'SIS', 'SERSIC', 'HERNQUIST']} @@ -144,6 +144,72 @@ def test_mass_scaling(self): assert kwargs_lens[1]['alpha_Rs'] == 0.1 assert kwargs_lens[2]['alpha_Rs'] == 0.3 + def test_general_scaling(self): + kwargs_model = {'lens_model_list': ['PJAFFE', 'PJAFFE', 'NFW', 'PJAFFE', 'NFW']} + # Scale Rs for two of the PJAFFEs, and sigma0 for a different set of PJAFFEs + # Scale alpha_Rs for the NFWs + kwargs_constraints = { + 'general_scaling': { + 'Rs': [1, False, False, 1, False], + 'sigma0': [False, 1, False, 1, False], + 'alpha_Rs': [False, False, 1, False, 1], + } + } + # PJAFFE: sigma0, Ra, Rs, center_x, center_y + # NFW: Rs, alpha_Rs, center_x, center_y + kwargs_fixed_lens = [ + {'Rs': 2.0, 'center_x': 1.0}, + {'sigma0': 2.0, 'Ra': 2.0, 'Rs': 3.0, 'center_y': 1.5}, + {'alpha_Rs': 0.1}, + {'Ra': 0.1, 'center_x': 0, 'center_y': 0}, + {'Rs': 3, 'center_x': -1, 'center_y': 3}, + ] + kwargs_fixed_cosmo = {} + param_class = Param(kwargs_model, kwargs_fixed_lens=kwargs_fixed_lens, kwargs_fixed_special=kwargs_fixed_cosmo + , **kwargs_constraints) + kwargs_lens = [{'sigma0': 3, 'Ra': 2, 'center_y': 5}, + {'center_x': 1.}, + {'Rs': 3, 'center_x': 0.0, 'center_y': -1.0}, + {'sigma0': 3, 'Rs': 1.5}, + {'alpha_Rs': 4}] + kwargs_source = [] + kwargs_lens_light = [] + kwargs_ps = [] + # Define the scaling and power for each parameter + kwargs_cosmo = { + 'Rs_scale_factor': [2.0], + 'Rs_scale_pow': [1.1], + 'sigma0_scale_factor': [3], + 'sigma0_scale_pow': [2.0], + 'alpha_Rs_scale_factor': [0.3], + 'alpha_Rs_scale_pow': [0.5], + } + args = param_class.kwargs2args(kwargs_lens, kwargs_source, kwargs_lens_light, kwargs_ps, kwargs_special=kwargs_cosmo) + num, names = param_class.num_param() + print(names) + print(args) + + kwargs_return = param_class.args2kwargs(args) + kwargs_lens = kwargs_return['kwargs_lens'] + print('kwargs_lens:', kwargs_lens) + npt.assert_almost_equal(kwargs_lens[0]['Rs'], 2.0 * 2.0**1.1) + npt.assert_almost_equal(kwargs_lens[0]['sigma0'], 3) + npt.assert_almost_equal(kwargs_lens[1]['Rs'], 3.0) + npt.assert_almost_equal(kwargs_lens[1]['sigma0'], 3.0 * 2.0**2.0) + npt.assert_almost_equal(kwargs_lens[2]['alpha_Rs'], 0.3 * 0.1**0.5) + npt.assert_almost_equal(kwargs_lens[3]['Rs'], 2.0 * 1.5**1.1) + npt.assert_almost_equal(kwargs_lens[3]['sigma0'], 3.0 * 3**2.0) + npt.assert_almost_equal(kwargs_lens[4]['alpha_Rs'], 0.3 * 4**0.5) + + kwargs_return = param_class.args2kwargs(args, bijective=True) + kwargs_lens = kwargs_return['kwargs_lens'] + npt.assert_almost_equal(kwargs_lens[0]['Rs'], 2.0) + npt.assert_almost_equal(kwargs_lens[1]['sigma0'], 2.0) + npt.assert_almost_equal(kwargs_lens[2]['alpha_Rs'], 0.1) + npt.assert_almost_equal(kwargs_lens[3]['Rs'], 1.5) + npt.assert_almost_equal(kwargs_lens[3]['sigma0'], 3) + npt.assert_almost_equal(kwargs_lens[4]['alpha_Rs'], 4) + def test_joint_lens_with_light(self): kwargs_model = {'lens_model_list': ['CHAMELEON'], 'lens_light_model_list': ['CHAMELEON']} i_light, k_lens = 0, 0 diff --git a/test/test_Sampling/test_special_param.py b/test/test_Sampling/test_special_param.py index 88edbfbee..68b8ded50 100644 --- a/test/test_Sampling/test_special_param.py +++ b/test/test_Sampling/test_special_param.py @@ -79,6 +79,18 @@ def test_source_grid_offsets(self): assert kwargs_new['delta_x_source_grid'] == kwargs_fixed['delta_x_source_grid'] assert kwargs_new['delta_y_source_grid'] == kwargs_fixed['delta_y_source_grid'] + def test_general_scaling(self): + kwargs_fixed = {} + param = SpecialParam(kwargs_fixed=kwargs_fixed, + general_scaling_params={'param': [False, 1, 1, False, 2]}) + args = param.set_params({'param_scale_factor': [1, 2], 'param_scale_pow': [3, 4]}) + assert len(args) == 4 + num_param, param_list = param.num_param() + assert num_param == 4 + kwargs_new, _ = param.get_params(args, i=0) + assert kwargs_new['param_scale_factor'] == [1, 2] + assert kwargs_new['param_scale_pow'] == [3, 4] + if __name__ == '__main__': pytest.main() diff --git a/test/test_SimulationAPI/test_data_api.py b/test/test_SimulationAPI/test_data_api.py index 9575acb9f..235fed430 100644 --- a/test/test_SimulationAPI/test_data_api.py +++ b/test/test_SimulationAPI/test_data_api.py @@ -30,6 +30,15 @@ def setup(self): kwargs_data = util.merge_dicts(kwargs_instrument, kwargs_observations) self.api_pixel = DataAPI(numpix=numpix, data_count_unit='ADU', **kwargs_data) + self.ra_at_xy_0 = 0.02 + self.dec_at_xy_0 = 0.02 + self.transform_pix2angle = [[-self.pixel_scale,0],[0,self.pixel_scale]] + kwargs_pixel_grid = {'ra_at_xy_0':self.ra_at_xy_0,'dec_at_xy_0':self.dec_at_xy_0, + 'transform_pix2angle':self.transform_pix2angle} + self.api_pixel_grid = DataAPI(numpix=numpix, + kwargs_pixel_grid=kwargs_pixel_grid, + data_count_unit='ADU',**self.kwargs_data) + def test_data_class(self): data_class = self.api.data_class assert data_class.pixel_width == self.pixel_scale @@ -40,6 +49,12 @@ def test_psf_class(self): psf_class = self.api_pixel.psf_class assert psf_class.psf_type == 'PIXEL' + def test_kwargs_data(self): + kwargs_data = self.api.kwargs_data + assert kwargs_data['ra_at_xy_0'] != self.ra_at_xy_0 + kwargs_data = self.api_pixel_grid.kwargs_data + assert kwargs_data['ra_at_xy_0'] == self.ra_at_xy_0 + class TestRaise(unittest.TestCase): @@ -72,3 +87,9 @@ def test_raise(self): with self.assertRaises(ValueError): data_api = DataAPI(numpix=numpix, data_count_unit='ADU', **kwargs_data) psf_class = data_api.psf_class + + kwargs_data['kernel_point_source'] = np.ones((3, 3)) + kwargs_pixel_grid = {'ra_at_xy_0':0.02,'dec_at_xy_0':0.02} + with self.assertRaises(ValueError): + data_api = DataAPI(numpix=numpix,kwargs_pixel_grid=kwargs_pixel_grid, + **kwargs_data) diff --git a/test/test_Util/test_analysis_util.py b/test/test_Util/test_analysis_util.py index d91c64e1e..c7e677dcf 100644 --- a/test/test_Util/test_analysis_util.py +++ b/test/test_Util/test_analysis_util.py @@ -47,8 +47,8 @@ def test_half_light_radius(self): assert r_half == -1 def test_bic_model(self): - bic=analysis_util.bic_model(0,np.e,1) - assert bic == 1 + bic=analysis_util.bic_model(0, np.e, 1) + npt.assert_almost_equal(bic, 1, decimal=8) def test_azimuthalAverage(self): num_pix = 101 diff --git a/test/test_Util/test_util.py b/test/test_Util/test_util.py index 9ca203e52..624e617fb 100644 --- a/test/test_Util/test_util.py +++ b/test/test_Util/test_util.py @@ -322,11 +322,22 @@ def test_min_square_dist(): assert dist[1] == 1 +def test_neighbor_select_fast(): + a = np.ones(100) + a[41] = 0 + x = np.linspace(0, 99, 100) + y = np.linspace(0, 99, 100) + x_mins, y_mins, values = util.local_minima_2d(a, x, y) + assert x_mins[0] == 41 + assert y_mins[0] == 41 + assert values[0] == 0 + + def test_neighborSelect(): a = np.ones(100) a[41] = 0 - x = np.linspace(0,99,100) - y = np.linspace(0,99,100) + x = np.linspace(0, 99, 100) + y = np.linspace(0, 99, 100) x_mins, y_mins, values = util.neighborSelect(a, x, y) assert x_mins[0] == 41 assert y_mins[0] == 41 diff --git a/test/test_Workflow/test_fitting_sequence.py b/test/test_Workflow/test_fitting_sequence.py index b5036e5b8..ed6f983ac 100644 --- a/test/test_Workflow/test_fitting_sequence.py +++ b/test/test_Workflow/test_fitting_sequence.py @@ -198,6 +198,7 @@ def test_fitting_sequence(self): assert kwargs_set['kwargs_ps'][0]['ra_source'] == 0.007 def test_zeus(self): + np.random.seed(42) # we make a very basic lens+source model to feed to check zeus can be run through fitting sequence # we don't use the kwargs defined in setup() as those are modified during the tests; using unique kwargs here is safer @@ -214,8 +215,6 @@ def test_zeus(self): data_class = ImageData(**kwargs_data) kwargs_psf_gaussian = {'psf_type': 'GAUSSIAN', 'fwhm': fwhm, 'pixel_size': deltaPix, 'truncation': 3} psf_gaussian = PSF(**kwargs_psf_gaussian) - kwargs_psf = {'psf_type': 'PIXEL', 'kernel_point_source': psf_gaussian.kernel_point_source, 'psf_error_map': np.zeros_like(psf_gaussian.kernel_point_source)} - psf_class = PSF(**kwargs_psf) # make a lens lens_model_list = ['EPL'] @@ -232,7 +231,7 @@ def test_zeus(self): kwargs_numerics = {'supersampling_factor': 1, 'supersampling_convolution': False} - imageModel = ImageModel(data_class, psf_class, lens_model_class, source_model_class, kwargs_numerics=kwargs_numerics) + imageModel = ImageModel(data_class, psf_gaussian, lens_model_class, source_model_class, kwargs_numerics=kwargs_numerics) image_sim = sim_util.simulate_simple(imageModel, kwargs_lens, kwargs_source) data_class.update_data(image_sim) @@ -260,12 +259,12 @@ def test_zeus(self): kwargs_constraints = {} - multi_band_list = [[kwargs_data, kwargs_psf, kwargs_numerics]] + multi_band_list = [[kwargs_data, kwargs_psf_gaussian, kwargs_numerics]] kwargs_data_joint = {'multi_band_list': multi_band_list, 'multi_band_type': 'multi-linear'} - kwargs_likelihood = {'source_marg': True} + kwargs_likelihood = {'source_marg': False} fittingSequence = FittingSequence(kwargs_data_joint, kwargs_model, kwargs_constraints, kwargs_likelihood, @@ -278,7 +277,6 @@ def test_zeus(self): chain_list = fittingSequence.fit_sequence(fitting_list) - def test_multinest(self): # Nested sampler tests # further decrease the parameter space for nested samplers to run faster @@ -321,8 +319,9 @@ def test_multinest(self): assert kwargs_out['kwargs_lens'] == 1 def test_dynesty(self): + np.random.seed(42) kwargs_params = copy.deepcopy(self.kwargs_params) - kwargs_params['lens_model'][0][0]['theta_E'] += 0.01 + kwargs_params['lens_model'][0][0]['theta_E'] += 0.2 fittingSequence = FittingSequence(self.kwargs_data_joint, self.kwargs_model, self.kwargs_constraints, self.kwargs_likelihood, kwargs_params) @@ -341,6 +340,7 @@ def test_dynesty(self): chain_list = fittingSequence.fit_sequence(fitting_list) def test_nautilus(self): + np.random.seed(42) kwargs_params = copy.deepcopy(self.kwargs_params) fittingSequence = FittingSequence(self.kwargs_data_joint, self.kwargs_model, self.kwargs_constraints, self.kwargs_likelihood, kwargs_params)