Merge pull request #536 from sibirrer/tracer_model

Tracer model more tests and in sync with main

.pre-commit-config.yaml

@@ -12,7 +12,7 @@ ci:
 
 repos:
   - repo: https://github.com/psf/black
-    rev: 23.7.0
+    rev: 23.9.1
     hooks:
       - id: black
         # It is recommended to specify the latest version of Python
@@ -21,7 +21,7 @@ repos:
         # https://pre-commit.com/#top_level-default_language_version
         language_version: python3
   - repo: https://github.com/psf/black
-    rev: 23.7.0
+    rev: 23.9.1
     hooks:
       - id: black-jupyter
         language_version: python3

.readthedocs.yml

@@ -5,6 +5,16 @@
 # Required
 version: 2
 
+# Set the OS, Python version and other tools you might need
+build:
+  os: ubuntu-22.04
+  tools:
+    python: "3.11"
+    # You can also specify other tool versions:
+    # nodejs: "20"
+    # rust: "1.70"
+    # golang: "1.20"
+
 # Build documentation in the docs/ directory with Sphinx
 sphinx:
   configuration: docs/conf.py
@@ -18,7 +28,5 @@ formats: all
 
 # Optionally set the version of Python and requirements required to build your docs
 python:
-  version: "3.8"
   install:
-    #- requirements: requirements.txt
     - requirements: requirements.txt

HISTORY.rst

@@ -348,4 +348,14 @@ History
 * cored truncated NFW profile
 * SkiNN interface
 * faster and more reliable Einstein radius computation
-* Cobaya interface
+* Cobaya interface
+
+1.11.4 (2023-09-27)
++++++++++++++++++++
+* Cobaya not required to run FittingSequence
+* Galkin with luminosity-weighted velocity dispersion calculation
+
+1.11.5 (2023-09-28)
++++++++++++++++++++
+* bug fix in findOverlap function
+* bug fix in luminosity-0weighted celocity dispersion calculation

lenstronomy/GalKin/numeric_kinematics.py

@@ -86,7 +86,7 @@ def lum_weighted_vel_disp(self, R, kwargs_mass, kwargs_light, kwargs_anisotropy)
         )
         # azimuthal averaging
         I_R_sigma2 = np.sum(I_R_sigma2_rad * r_rad)
-        I_R = np.sum(I_R_rad)
+        I_R = np.sum(I_R_rad * r_rad)
         # return in units km/s
         sigma_s2 = I_R_sigma2 / I_R
         sigma = np.sqrt(sigma_s2) / 1000.0

lenstronomy/LensModel/Profiles/base_profile.py

@@ -4,7 +4,25 @@
 class LensProfileBase(object):
     """This class acts as the base class of all lens model functions and indicates raise
     statements and default outputs if these functions are not defined in the specific
-    lens model class."""
+    lens model class.
+
+    To implement a new lens model profile you should:
+    1. make a new python file in this folder
+    2. create a class inheriting this class; YourModel(LensProfileBase)
+    3. write new definitions following the same input and output conventions as this base class
+        function(x, y, <other parameters>)
+        derivatives(x, y, <other parameters>)
+        hessian(x, y, <other parameters>)
+    4. set the variables for sampling the new profile
+        param_names = ["param1", "param2", ...]
+        lower_limit_default = {"param1": value, "param2: value, ...}
+        upper_limit_default = {"param1": value, "param2: value, ...}
+    5. give the new profile a meaningful name and add it in the LensModel.profile_list_base class
+    6. write test functions in the test/test_LensModel/test_Profiles folder with a new file with test_<profile name>.py
+
+    With that, you should be good to go and import and use it for any purpose.
+    Further definitions in the class are optional and only used for certain applications (such as kinematics)
+    """
 
     def __init__(self, *args, **kwargs):
         self._static = False

lenstronomy/LensModel/Profiles/nfw.py

@@ -250,7 +250,7 @@ def nfwAlpha(self, R, Rs, rho0, ax_x, ax_y):
         R = np.maximum(R, 0.00000001)
         x = R / Rs
         gx = self.g_(x)
-        a = 4 * rho0 * Rs * R * gx / x**2 / R
+        a = 4 * rho0 * Rs * gx / x**2
         return a * ax_x, a * ax_y
 
     def nfwGamma(self, R, Rs, rho0, ax_x, ax_y):

lenstronomy/LensModel/Profiles/radial_interpolated.py

@@ -0,0 +1,199 @@
+import numpy as np
+from lenstronomy.LensModel.Profiles.base_profile import LensProfileBase
+from lenstronomy.Util import param_util
+from scipy.interpolate import interp1d
+from scipy import integrate
+
+__all__ = ["RadialInterpolate"]
+
+
+class RadialInterpolate(LensProfileBase):
+    """Radially interpolated profile with azimuthal symmetry."""
+
+    param_names = [
+        "r_bin",
+        "kappa_r",
+        "center_x",
+        "center_y",
+    ]
+    lower_limit_default = {}
+    upper_limit_default = {}
+
+    def function(
+        self,
+        x,
+        y,
+        r_bin=None,
+        kappa_r=None,
+        center_x=0,
+        center_y=0,
+    ):
+        """Lensing potential (only needed for specific calculations, such as time
+        delays)
+
+        :param x: x-coordinate (angular position), float or numpy array
+        :param y: y-coordinate (angular position), float or numpy array
+        :param r_bin: radial bins for which convergence values are provided
+        :type r_bin: array
+        :param kappa_r: convergence values corresponding to the r_bin radii
+        :type kappa_r: array of same size as r_bin
+        :param center_x: x-position of center of radial density profile
+        :type center_x: float
+        :param center_y: y-position of center of radial density profile
+        :type center_y: float
+        :return: lensing potential
+        """
+        r, phi = param_util.cart2polar(x, y, center_x=center_x, center_y=center_y)
+        # -\int alpha(r) dr  from 0 to r
+        pot = self._potential_r(r, r_bin, kappa_r)
+        return pot
+
+    def derivatives(
+        self,
+        x,
+        y,
+        r_bin=None,
+        kappa_r=None,
+        center_x=0,
+        center_y=0,
+    ):
+        """Returns df/dx and df/dy of the function.
+
+        :param x: x-coordinate (angular position), float or numpy array
+        :param y: y-coordinate (angular position), float or numpy array
+        :param r_bin: radial bins for which convergence values are provided
+        :type r_bin: array
+        :param kappa_r: convergence values corresponding to the r_bin radii
+        :type kappa_r: array of same size as r_bin
+        :param center_x: x-position of center of radial density profile
+        :type center_x: float
+        :param center_y: y-position of center of radial density profile
+        :type center_y: float
+        :return: f_x, f_y at interpolated positions (x, y)
+        """
+        x_ = x - center_x
+        y_ = y - center_y
+        r = np.maximum(np.sqrt(x_**2 + y_**2), 10 ** (-10))
+        alpha = self.alpha(r, r_bin, kappa_r)
+        return alpha * x_ / r, alpha * y_ / r
+
+    def hessian(
+        self,
+        x,
+        y,
+        r_bin=None,
+        kappa_r=None,
+        center_x=0,
+        center_y=0,
+    ):
+        """Returns Hessian matrix of function d^2f/dx^2, d^2/dxdy, d^2/dydx, d^f/dy^2.
+
+        :param x: x-coordinate (angular position), float or numpy array
+        :param y: y-coordinate (angular position), float or numpy array
+        :param r_bin: radial bins for which convergence values are provided
+        :type r_bin: array
+        :param kappa_r: convergence values corresponding to the r_bin radii
+        :type kappa_r: array of same size as r_bin
+        :param center_x: x-position of center of radial density profile
+        :type center_x: float
+        :param center_y: y-position of center of radial density profile
+        :type center_y: float
+        :return: f_xx, f_xy, f_yx, f_yy at interpolated positions (x, y)
+        """
+        r, phi = param_util.cart2polar(x, y, center_x=center_x, center_y=center_y)
+        r = np.maximum(r, 10 ** (-10))
+        kappa = self._kappa_r_interp(r, r_bin, kappa_r)
+        # shear in the spherical case is the average convergence enclosed minus the convergence at the radius
+        # source: Kaiser 1995
+        gamma = 1 * (self._mass_enclosed(r, r_bin, kappa) / (np.pi * r**2) - kappa)
+        # turn in to vector for gamma and cartesian derivatives
+        gamma1 = -np.cos(2.0 * phi) * gamma
+        gamma2 = -np.sin(2.0 * phi) * gamma
+        f_xx = kappa + gamma1
+        f_yy = kappa - gamma1
+        f_xy = gamma2
+
+        return f_xx, f_xy, f_xy, f_yy
+
+    def alpha(self, r, r_bin, kappa_r):
+        """Radial deflection angle m(<r) / r / pi.
+
+        :param r: radius from center
+        :param r_bin: radial bins for which convergence values are provided
+        :type r_bin: array
+        :param kappa_r: convergence values corresponding to the r_bin radii
+        :type kappa_r: array of same size as r_bin
+        :return: radial deflection angle
+        """
+        r_ = np.maximum(r, 10 ** (-10))
+        return self._mass_enclosed(r, r_bin, kappa_r) / r_ / np.pi
+
+    def _kappa_r_interp(self, r, r_bin, kappa_r):
+        """Calls interpolated kappa(r)
+
+        :param r: radius
+        :param r_bin: radial bins for which convergence values are provided
+        :type r_bin: numpy array
+        :param kappa_r: convergence values corresponding to the r_bin radii
+        :type kappa_r: array of same size as r_bin
+        :return: kappa(r)
+        """
+        if not hasattr(self, "_interp_kappa"):
+            self._interp_kappa = interp1d(
+                r_bin, kappa_r, fill_value=(kappa_r[0], kappa_r[-1]), bounds_error=False
+            )
+        return self._interp_kappa(r)
+
+    def _mass_enclosed(self, r, r_bin, kappa_r):
+        """Convergence enclosed a radius.
+
+        :param r: radius
+        :param r_bin: radial bins for which convergence values are provided
+        :type r_bin: numpy array
+        :param kappa_r: convergence values corresponding to the r_bin radii
+        :type kappa_r: array of same size as r_bin
+        :return: integrated convergence within radius <r
+        """
+        if not hasattr(self, "_interp_m_enclosed"):
+
+            def _integrand(x):
+                return x * 2 * np.pi * self._kappa_r_interp(x, r_bin, kappa_r)
+
+            m_slice_list = []
+            r_min = 0
+            for r_ in r_bin:
+                m_slice, _ = integrate.quad(_integrand, r_min, r_)
+                m_slice_list.append(m_slice)
+                r_min = r_
+            m_r = np.cumsum(m_slice_list)
+            self._interp_m_enclosed = interp1d(
+                r_bin, m_r, fill_value=(0, m_r[-1]), bounds_error=False
+            )
+        return self._interp_m_enclosed(r)
+
+    def _potential_r(self, r, r_bin, kappa_r):
+        """Convergence enclosed a radius.
+
+        :param r: radius
+        :param r_bin: radial bins for which convergence values are provided
+        :type r_bin: numpy array
+        :param kappa_r: convergence values corresponding to the r_bin radii
+        :type kappa_r: array of same size as r_bin
+        :return: integrated convergence within radius <r
+        """
+        if not hasattr(self, "_interp_potential"):
+
+            def _integrand(x):
+                return self.alpha(x, r_bin, kappa_r)
+
+            pot_slice_list = []
+            r_min = 0
+            for r_ in r_bin:
+                pot_slice, _ = integrate.quad(_integrand, r_min, r_)
+                pot_slice_list.append(pot_slice)
+                r_min = r_
+            pot_r = np.cumsum(pot_slice_list)
+            self._interp_potential = interp1d(
+                r_bin, pot_r, fill_value=(0, pot_r[-1]), bounds_error=False
+            )
+        return self._interp_potential(r)

lenstronomy/LensModel/profile_list_base.py

@@ -45,7 +45,7 @@
     "HESSIAN",
     "INTERPOL",
     "INTERPOL_SCALED",
-    "LOS",
+    "RADIAL_INTERPOL" "LOS",
     "LOS_MINIMAL",
     "MULTIPOLE",
     "MULTI_GAUSSIAN_KAPPA",
@@ -156,7 +156,7 @@ def _load_model_instances(
                 "CTNFW_GAUSS_DEC",
                 "INTERPOL",
                 "INTERPOL_SCALED",
-                "NIE",
+                "RADIAL_INTERPOL" "NIE",
                 "NIE_SIMPLE",
             ]:
                 lensmodel_class = self._import_class(
@@ -447,6 +447,14 @@ def _import_class(
             from lenstronomy.LensModel.Profiles.interpol import InterpolScaled
 
             return InterpolScaled(**kwargs_interp)
+
+        elif lens_type == "RADIAL_INTERPOL":
+            from lenstronomy.LensModel.Profiles.radial_interpolated import (
+                RadialInterpolate,
+            )
+
+            return RadialInterpolate()
+
         elif lens_type == "SHAPELETS_POLAR":
             from lenstronomy.LensModel.Profiles.shapelet_pot_polar import PolarShapelets
 

lenstronomy/Util/image_util.py

@@ -180,8 +180,8 @@ def findOverlap(x_mins, y_mins, min_distance):
             pass
         else:
             for j in range(0, i):
-                if abs(
-                    x_mins[i] - x_mins[j] < min_distance
+                if (
+                    abs(x_mins[i] - x_mins[j]) < min_distance
                     and abs(y_mins[i] - y_mins[j]) < min_distance
                 ):
                     idex.append(i)

lenstronomy/Util/util.py

@@ -233,6 +233,32 @@ def make_grid_transformed(numPix, Mpix2Angle):
     return ra_grid, dec_grid
 
 
+def centered_coordinate_system(num_pix, transform_pix2angle):
+    """Returns dictionary for Coordinate Grid such that (0,0) is centered with given
+    input orientation coordinate transformation matrix.
+
+    :param num_pix: number of pixels
+    :type num_pix: int
+    :param transform_pix2angle: transformation matrix (2x2) of pixels into coordinate
+        displacements
+    :return: dict with ra_at_xy_0, dec_at_xy_0, transfrom_pix2angle
+    """
+    pix_center = (num_pix - 1) / 2
+    ra_center = (
+        pix_center * transform_pix2angle[0, 0] + pix_center * transform_pix2angle[0, 1]
+    )
+    dec_center = (
+        pix_center * transform_pix2angle[1, 0] + pix_center * transform_pix2angle[1, 1]
+    )
+
+    kwargs_grid = {
+        "ra_at_xy_0": -ra_center,
+        "dec_at_xy_0": -dec_center,
+        "transform_pix2angle": transform_pix2angle,
+    }
+    return kwargs_grid
+
+
 @export
 def make_grid_with_coordtransform(
     numPix,

lenstronomy/__init__.py

@@ -1,6 +1,6 @@
 __author__ = "lenstronomy developers"
 __email__ = "lenstronomy-dev@googlegroups.com"
-__version__ = "1.11.3"
+__version__ = "1.11.5"
 __credits__ = "ETH Zurich, UCLA, Stanford, Stony Brook"
 
 from .Util.package_util import short, laconic