Repository: lenstronomy/lenstronomy (language: python)

Issue:
PyCosmos dependence
remove PyCosmo dependency to let test functions work


Discussion / hints:
pycosmo dependence removed and replaced with astropy.cosmology

Fix patch:
diff --git a/lenstronomy/Data/kinematic_bin_2D.py b/lenstronomy/Data/kinematic_bin_2D.py
new file mode 100644
index 000000000..fe25b481a
--- /dev/null
+++ b/lenstronomy/Data/kinematic_bin_2D.py
@@ -0,0 +1,77 @@
+import numpy as np
+import lenstronomy.Util.util as util
+from lenstronomy.Data.pixel_grid import PixelGrid
+
+__all__ = ['KinBin']
+
+
+class KinBin(object):
+    """
+    Class that summarizes the binned kinematic data.
+
+    The KinBin() class is initialized with :
+     - The information about the bins (bin values, and bin covariances, which pixels belong to which bin):
+    bin_data, bin_cov, bin_mask.
+     - The information about the associated intial shape of the unbinned kinematic map: bin_mask gives the index of
+    corresponding bin for each pixel), and ra_at_xy_0,dec_at_xy_0,transform_pix2angle,ra_shift,dec_shift are the usual
+    PixelGrid characteritics.
+
+    """
+
+    def __init__(self, bin_data, bin_cov, bin_mask, ra_at_xy_0, dec_at_xy_0,
+                 transform_pix2angle, psf_class, ra_shift=0, dec_shift=0):
+        """
+        :param bin_data: list, kinematic value of each bin, ordered by bin index.
+        :param bin_cov: 2D array (nbins x nbins), vrms covariance matrix associated to each bin, ordered by bin index
+        :param bin_mask: 2D array, mapping from the unbinned image to the binned one, each pixel value is the
+         corresponding bin index.
+        :param ra_at_xy_0: float, ra coordinate at pixel (0,0) (unbinned image)
+        :param dec_at_xy_0: float, dec coordinate at pixel (0,0) (unbinned image)
+        :param transform_pix2angle: 2x2 array, mapping of pixel (unbinned image) to coordinate
+        :param psf_class: PSF class
+        :param ra_shift:  float, RA shift of pixel grid
+        :param dec_shift: float, DEC shift of pixel grid
+
+        """
+        self.PSF = psf_class
+        nx, ny = np.shape(bin_mask)
+        self._nx = nx
+        self._ny = ny
+        self.PixelGrid = PixelGrid(nx, ny, transform_pix2angle, ra_at_xy_0 + ra_shift, dec_at_xy_0 + dec_shift)
+
+        self.data = bin_data
+        self.covariance = bin_cov
+        self.bin_mask = bin_mask
+        self._pix2a = transform_pix2angle
+        self._ra_at_xy_0 = ra_at_xy_0
+        self._dec_at_xy_0 = dec_at_xy_0
+
+    @staticmethod
+    def binned_image(data, bin_mask):
+        """
+        Creates the binned image of the data
+
+        :param data: data value in each bin
+        :param bin_mask: mask indicating which pixels belong to which bin
+        """
+        binned_image = np.zeros_like(bin_mask)
+        for idx, value in enumerate(data):
+            binned_image[bin_mask == idx] = value
+        return binned_image
+
+    def kin_bin2kwargs(self):
+        """
+        Creates the kwargs needed for the 2D kinematic likelihood
+        """
+        kwargs = {'image': self.binned_image(self.data, self.bin_mask),
+                  'deltaPix': self.PixelGrid.pixel_width,
+                  'transform_pix2angle': self._pix2a, 'ra_at_xy0': self._ra_at_xy_0,
+                  'dec_at_xy0': self._dec_at_xy_0}
+        return kwargs
+
+    def kin_grid(self):
+        """
+        Creates a pixel grid that satisfy the kinematics coordinates system
+        """
+        x_grid, y_grid = self.PixelGrid.pixel_coordinates
+        return x_grid, y_grid
diff --git a/lenstronomy/LensModel/Profiles/epl.py b/lenstronomy/LensModel/Profiles/epl.py
index d6ba158c5..5c793c18b 100644
--- a/lenstronomy/LensModel/Profiles/epl.py
+++ b/lenstronomy/LensModel/Profiles/epl.py
@@ -7,8 +7,7 @@
 from lenstronomy.LensModel.Profiles.spp import SPP
 from scipy.special import hyp2f1
 
-
-__all__ = ['EPL', 'EPLMajorAxis']
+__all__ = ['EPL', 'EPLMajorAxis', 'EPLQPhi']
 
 
 class EPL(LensProfileBase):
@@ -188,8 +187,8 @@ def hessian(self, x, y, theta_E, gamma, e1, e2, center_x=0, center_y=0):
         # evaluate
         f__xx, f__xy, f__yx, f__yy = self.epl_major_axis.hessian(x__, y__, b, t, q)
         # rotate back
-        kappa = 1./2 * (f__xx + f__yy)
-        gamma1__ = 1./2 * (f__xx - f__yy)
+        kappa = 1. / 2 * (f__xx + f__yy)
+        gamma1__ = 1. / 2 * (f__xx - f__yy)
         gamma2__ = f__xy
         gamma1 = np.cos(2 * phi_G) * gamma1__ - np.sin(2 * phi_G) * gamma2__
         gamma2 = +np.sin(2 * phi_G) * gamma1__ + np.cos(2 * phi_G) * gamma2__
@@ -241,7 +240,6 @@ class EPLMajorAxis(LensProfileBase):
     param_names = ['b', 't', 'q', 'center_x', 'center_y']
 
     def __init__(self):
-
         super(EPLMajorAxis, self).__init__()
 
     def function(self, x, y, b, t, q):
@@ -259,7 +257,7 @@ def function(self, x, y, b, t, q):
         alpha_x, alpha_y = self.derivatives(x, y, b, t, q)
 
         # deflection potential, eq. (15)
-        psi = (x*alpha_x + y*alpha_y)/(2 - t)
+        psi = (x * alpha_x + y * alpha_y) / (2 - t)
 
         return psi
 
@@ -276,20 +274,20 @@ def derivatives(self, x, y, b, t, q):
         """
         # elliptical radius, eq. (5)
         Z = np.empty(np.shape(x), dtype=complex)
-        Z.real = q*x
+        Z.real = q * x
         Z.imag = y
         R = np.abs(Z)
         R = np.maximum(R, 0.000000001)
 
         # angular dependency with extra factor of R, eq. (23)
-        R_omega = Z*hyp2f1(1, t/2, 2-t/2, -(1-q)/(1+q)*(Z/Z.conj()))
+        R_omega = Z * hyp2f1(1, t / 2, 2 - t / 2, -(1 - q) / (1 + q) * (Z / Z.conj()))
 
         # deflection, eq. (22)
-        alpha = 2/(1+q)*(b/R)**t*R_omega
+        alpha = 2 / (1 + q) * (b / R) ** t * R_omega
 
         # return real and imaginary part
-        alpha_real = np.nan_to_num(alpha.real, posinf=10**10, neginf=-10**10)
-        alpha_imag = np.nan_to_num(alpha.imag, posinf=10**10, neginf=-10**10)
+        alpha_real = np.nan_to_num(alpha.real, posinf=10 ** 10, neginf=-10 ** 10)
+        alpha_imag = np.nan_to_num(alpha.imag, posinf=10 ** 10, neginf=-10 ** 10)
 
         return alpha_real, alpha_imag
 
@@ -304,25 +302,25 @@ def hessian(self, x, y, b, t, q):
         :param q: axis ratio
         :return: f_xx, f_yy, f_xy
         """
-        R = np.hypot(q*x, y)
+        R = np.hypot(q * x, y)
         R = np.maximum(R, 0.00000001)
         r = np.hypot(x, y)
 
-        cos, sin = x/r, y/r
-        cos2, sin2 = cos*cos*2 - 1, sin*cos*2
+        cos, sin = x / r, y / r
+        cos2, sin2 = cos * cos * 2 - 1, sin * cos * 2
 
         # convergence, eq. (2)
-        kappa = (2 - t)/2*(b/R)**t
-        kappa = np.nan_to_num(kappa, posinf=10**10, neginf=-10**10)
+        kappa = (2 - t) / 2 * (b / R) ** t
+        kappa = np.nan_to_num(kappa, posinf=10 ** 10, neginf=-10 ** 10)
 
         # deflection via method
         alpha_x, alpha_y = self.derivatives(x, y, b, t, q)
 
         # shear, eq. (17), corrected version from arXiv/corrigendum
-        gamma_1 = (1-t)*(alpha_x*cos - alpha_y*sin)/r - kappa*cos2
-        gamma_2 = (1-t)*(alpha_y*cos + alpha_x*sin)/r - kappa*sin2
-        gamma_1 = np.nan_to_num(gamma_1, posinf=10**10, neginf=-10**10)
-        gamma_2 = np.nan_to_num(gamma_2, posinf=10**10, neginf=-10**10)
+        gamma_1 = (1 - t) * (alpha_x * cos - alpha_y * sin) / r - kappa * cos2
+        gamma_2 = (1 - t) * (alpha_y * cos + alpha_x * sin) / r - kappa * sin2
+        gamma_1 = np.nan_to_num(gamma_1, posinf=10 ** 10, neginf=-10 ** 10)
+        gamma_2 = np.nan_to_num(gamma_2, posinf=10 ** 10, neginf=-10 ** 10)
 
         # second derivatives from convergence and shear
         f_xx = kappa + gamma_1
@@ -330,3 +328,91 @@ def hessian(self, x, y, b, t, q):
         f_xy = gamma_2
 
         return f_xx, f_xy, f_xy, f_yy
+
+
+class EPLQPhi(LensProfileBase):
+    """
+    class to model a EPL sampling over q and phi instead of e1 and e2.
+
+    """
+    param_names = ['theta_E', 'gamma', 'q', 'phi', 'center_x', 'center_y']
+    lower_limit_default = {'theta_E': 0, 'gamma': 1.5, 'q': 0, 'phi': -np.pi, 'center_x': -100, 'center_y': -100}
+    upper_limit_default = {'theta_E': 100, 'gamma': 2.5, 'q': 1, 'phi': np.pi, 'center_x': 100, 'center_y': 100}
+
+    def __init__(self):
+        self._EPL = EPL()
+        super(EPLQPhi, self).__init__()
+
+    def function(self, x, y, theta_E, gamma, q, phi, center_x=0, center_y=0):
+        """
+
+        :param x: x-coordinate in image plane
+        :param y: y-coordinate in image plane
+        :param theta_E: Einstein radius
+        :param gamma: power law slope
+        :param q: axis ratio
+        :param phi: position angle
+        :param center_x: profile center
+        :param center_y: profile center
+        :return: lensing potential
+        """
+        e1, e2 = param_util.phi_q2_ellipticity(phi, q)
+        return self._EPL.function(x, y, theta_E, gamma, e1, e2, center_x, center_y)
+
+    def derivatives(self, x, y, theta_E, gamma, q, phi, center_x=0, center_y=0):
+        """
+
+        :param x: x-coordinate in image plane
+        :param y: y-coordinate in image plane
+        :param theta_E: Einstein radius
+        :param gamma: power law slope
+        :param q: axis ratio
+        :param phi: position angle
+        :param center_x: profile center
+        :param center_y: profile center
+        :return: alpha_x, alpha_y
+        """
+        e1, e2 = param_util.phi_q2_ellipticity(phi, q)
+        return self._EPL.derivatives(x, y, theta_E, gamma, e1, e2, center_x, center_y)
+
+    def hessian(self, x, y, theta_E, gamma, q, phi, center_x=0, center_y=0):
+        """
+
+        :param x: x-coordinate in image plane
+        :param y: y-coordinate in image plane
+        :param theta_E: Einstein radius
+        :param gamma: power law slope
+        :param q: axis ratio
+        :param phi: position angle
+        :param center_x: profile center
+        :param center_y: profile center
+        :return: f_xx, f_xy, f_yx, f_yy
+        """
+        e1, e2 = param_util.phi_q2_ellipticity(phi, q)
+        return self._EPL.hessian(x, y, theta_E, gamma, e1, e2, center_x, center_y)
+
+    def mass_3d_lens(self, r, theta_E, gamma, q=None, phi=None):
+        """
+        computes the spherical power-law mass enclosed (with SPP routine)
+        :param r: radius within the mass is computed
+        :param theta_E: Einstein radius
+        :param gamma: power-law slope
+        :param q: axis ratio (not used)
+        :param phi: position angle (not used)
+        :return: mass enclosed a 3D radius r
+        """
+        return self._EPL.mass_3d_lens(r, theta_E, gamma)
+
+    def density_lens(self, r, theta_E, gamma, q=None, phi=None):
+        """
+        computes the density at 3d radius r given lens model parameterization.
+        The integral in the LOS projection of this quantity results in the convergence quantity.
+
+        :param r: radius within the mass is computed
+        :param theta_E: Einstein radius
+        :param gamma: power-law slope
+        :param q: axis ratio (not used)
+        :param phi: position angle (not used)
+        :return: mass enclosed a 3D radius r
+        """
+        return self._EPL.density_lens(r, theta_E, gamma)
diff --git a/lenstronomy/LensModel/Profiles/pemd.py b/lenstronomy/LensModel/Profiles/pemd.py
index da3962cf4..5f363db6e 100644
--- a/lenstronomy/LensModel/Profiles/pemd.py
+++ b/lenstronomy/LensModel/Profiles/pemd.py
@@ -3,13 +3,15 @@
 from lenstronomy.LensModel.Profiles.spp import SPP
 from lenstronomy.LensModel.Profiles.spemd import SPEMD
 from lenstronomy.LensModel.Profiles.base_profile import LensProfileBase
+import lenstronomy.Util.param_util as param_util
+import numpy as np
 
 __all__ = ['PEMD']
 
 
 class PEMD(LensProfileBase):
     """
-    class for power law ellipse mass density profile.
+    class for power law ellipse mass density profile (PEMD).
     This class effectively calls the class SPEMD_SMOOTH with a fixed and very small central smoothing scale
     to perform the numerical integral using the FASTELL code by Renan Barkana.
     An alternative implementation of the same model using pure python with analytical functions is probided as 'EPL'
@@ -123,4 +125,4 @@ def density_lens(self, r, theta_E, gamma, e1=None, e2=None):
         :param e2: eccentricity component (not used)
         :return: mass enclosed a 3D radius r
         """
-        return self.spp.density_lens(r, theta_E, gamma)
+        return self.spp.density_lens(r, theta_E, gamma)
\ No newline at end of file
diff --git a/lenstronomy/LensModel/profile_list_base.py b/lenstronomy/LensModel/profile_list_base.py
index ae7969ba5..fa96be82d 100644
--- a/lenstronomy/LensModel/profile_list_base.py
+++ b/lenstronomy/LensModel/profile_list_base.py
@@ -5,7 +5,7 @@
 
 _SUPPORTED_MODELS = ['SHIFT', 'NIE_POTENTIAL', 'CONST_MAG', 'SHEAR', 'SHEAR_GAMMA_PSI', 'SHEAR_REDUCED', 'CONVERGENCE', 'FLEXION',
                      'FLEXIONFG', 'POINT_MASS', 'SIS', 'SIS_TRUNCATED', 'SIE', 'SPP', 'NIE', 'NIE_SIMPLE', 'CHAMELEON',
-                     'DOUBLE_CHAMELEON', 'TRIPLE_CHAMELEON', 'SPEP', 'PEMD', 'SPEMD', 'EPL', 'EPL_NUMBA',
+                     'DOUBLE_CHAMELEON', 'TRIPLE_CHAMELEON', 'SPEP', 'PEMD', 'SPEMD', 'EPL', 'EPL_Q_PHI', 'EPL_NUMBA',
                      'EPL_BOXYDISKY', 'SPL_CORE',
                      'NFW', 'NFW_ELLIPSE', 'NFW_ELLIPSE_GAUSS_DEC', 'NFW_ELLIPSE_CSE', 'TNFW', 'TNFW_ELLIPSE',
                      'CNFW', 'CNFW_ELLIPSE', 'CTNFW_GAUSS_DEC', 'NFW_MC', 'SERSIC',
@@ -159,6 +159,9 @@ def _import_class(lens_type, custom_class, kwargs_interp, kwargs_synthesis, z_le
         elif lens_type == 'EPL':
             from lenstronomy.LensModel.Profiles.epl import EPL
             return EPL()
+        elif lens_type == 'EPL_Q_PHI':
+            from lenstronomy.LensModel.Profiles.epl import EPLQPhi
+            return EPLQPhi()
         elif lens_type == 'EPL_NUMBA':
             from lenstronomy.LensModel.Profiles.epl_numba import EPL_numba
             return EPL_numba()
diff --git a/lenstronomy/LightModel/Profiles/sersic.py b/lenstronomy/LightModel/Profiles/sersic.py
index ab840274d..6fee5e1c8 100644
--- a/lenstronomy/LightModel/Profiles/sersic.py
+++ b/lenstronomy/LightModel/Profiles/sersic.py
@@ -4,7 +4,8 @@
 
 import numpy as np
 from lenstronomy.LensModel.Profiles.sersic_utils import SersicUtil
-
+import lenstronomy.Util.param_util as param_util
+import numpy as np
 from lenstronomy.Util.package_util import exporter
 export, __all__ = exporter()
 
@@ -87,6 +88,40 @@ def function(self, x, y, amp, R_sersic, n_sersic, e1, e2, center_x=0, center_y=0
         result = self._r_sersic(R, R_sersic, n_sersic, max_R_frac)
         return amp * result
 
+@export
+class SersicElliptic_qPhi(SersicUtil):
+    """
+    this class is the same as SersicElliptic except sampling over q and phi instead of e1 and e2
+
+    """
+    param_names = ['amp', 'R_sersic', 'n_sersic', 'q', 'phi', 'center_x', 'center_y']
+    lower_limit_default = {'amp': 0, 'R_sersic': 0, 'n_sersic': 0.5, 'q': 0, 'phi': -np.pi, 'center_x': -100,
+                           'center_y': -100}
+    upper_limit_default = {'amp': 100, 'R_sersic': 100, 'n_sersic': 8, 'q': 1.0, 'phi': np.pi, 'center_x': 100,
+                           'center_y': 100}
+
+    def __init__(self, *args, **kwargs):
+        self._sersic_e1e2 = SersicElliptic(*args, **kwargs)
+
+    def function(self, x, y, amp, R_sersic, n_sersic, q, phi, center_x=0, center_y=0, max_R_frac=100.0):
+        """
+
+        :param x:
+        :param y:
+        :param amp: surface brightness/amplitude value at the half light radius
+        :param R_sersic: half light radius (either semi-major axis or product average of semi-major and semi-minor axis)
+        :param n_sersic: Sersic index
+        :param q: axis ratio
+        :param phi: position angle (radians)
+        :param center_x: center in x-coordinate
+        :param center_y: center in y-coordinate
+        :param max_R_frac: maximum window outside of which the mass is zeroed, in units of R_sersic (float)
+        :return: Sersic profile value at (x, y)
+        """
+
+        e1, e2 = param_util.phi_q2_ellipticity(phi, q)
+        return self._sersic_e1e2.function(x, y, amp, R_sersic, n_sersic, e1, e2, center_x, center_y, max_R_frac)
+
 
 @export
 class CoreSersic(SersicUtil):
diff --git a/lenstronomy/LightModel/light_model_base.py b/lenstronomy/LightModel/light_model_base.py
index 79fd958b7..3b265c733 100644
--- a/lenstronomy/LightModel/light_model_base.py
+++ b/lenstronomy/LightModel/light_model_base.py
@@ -12,7 +12,7 @@
 
 
 _MODELS_SUPPORTED = ['GAUSSIAN', 'GAUSSIAN_ELLIPSE', 'ELLIPSOID', 'MULTI_GAUSSIAN', 'MULTI_GAUSSIAN_ELLIPSE',
-                     'SERSIC', 'SERSIC_ELLIPSE', 'CORE_SERSIC', 'SHAPELETS', 'SHAPELETS_POLAR', 'SHAPELETS_POLAR_EXP',
+                     'SERSIC', 'SERSIC_ELLIPSE', 'SERSIC_ELLIPSE_Q_PHI', 'CORE_SERSIC', 'SHAPELETS', 'SHAPELETS_POLAR', 'SHAPELETS_POLAR_EXP',
                      'SHAPELETS_ELLIPSE',
                      'HERNQUIST', 'HERNQUIST_ELLIPSE', 'PJAFFE', 'PJAFFE_ELLIPSE', 'UNIFORM', 'POWER_LAW', 'NIE',
                      'CHAMELEON', 'DOUBLE_CHAMELEON', 'TRIPLE_CHAMELEON', 'INTERPOL', 'SLIT_STARLETS',
@@ -58,6 +58,9 @@ def __init__(self, light_model_list, smoothing=0.001, sersic_major_axis=None):
             elif profile_type == 'SERSIC_ELLIPSE':
                 from lenstronomy.LightModel.Profiles.sersic import SersicElliptic
                 self.func_list.append(SersicElliptic(smoothing=smoothing, sersic_major_axis=sersic_major_axis))
+            elif profile_type == 'SERSIC_ELLIPSE_Q_PHI':
+                from lenstronomy.LightModel.Profiles.sersic import SersicElliptic_qPhi
+                self.func_list.append(SersicElliptic_qPhi(smoothing=smoothing, sersic_major_axis=sersic_major_axis))
             elif profile_type == 'CORE_SERSIC':
                 from lenstronomy.LightModel.Profiles.sersic import CoreSersic
                 self.func_list.append(CoreSersic(smoothing=smoothing, sersic_major_axis=sersic_major_axis))
diff --git a/lenstronomy/LightModel/linear_basis.py b/lenstronomy/LightModel/linear_basis.py
index 720878b43..d54fb1229 100644
--- a/lenstronomy/LightModel/linear_basis.py
+++ b/lenstronomy/LightModel/linear_basis.py
@@ -61,7 +61,8 @@ def functions_split(self, x, y, kwargs_list, k=None):
         bool_list = self._bool_list(k=k)
         for i, model in enumerate(self.profile_type_list):
             if bool_list[i] is True:
-                if model in ['SERSIC', 'SERSIC_ELLIPSE', 'CORE_SERSIC', 'HERNQUIST', 'HERNQUIST_ELLIPSE', 'PJAFFE',
+                if model in ['SERSIC', 'SERSIC_ELLIPSE', 'SERSIC_ELLIPSE_Q_PHI', 'CORE_SERSIC',
+                             'HERNQUIST', 'HERNQUIST_ELLIPSE', 'PJAFFE',
                              'PJAFFE_ELLIPSE', 'GAUSSIAN', 'GAUSSIAN_ELLIPSE', 'POWER_LAW', 'NIE', 'CHAMELEON',
                              'DOUBLE_CHAMELEON', 'TRIPLE_CHAMELEON', 'UNIFORM', 'INTERPOL', 'ELLIPSOID']:
                     kwargs_new = kwargs_list[i].copy()
@@ -150,7 +151,8 @@ def update_linear(self, param, i, kwargs_list):
         :return: kwargs list with over-written or added 'amp' parameters according to the coefficients in param
         """
         for k, model in enumerate(self.profile_type_list):
-            if model in ['SERSIC', 'SERSIC_ELLIPSE', 'CORE_SERSIC', 'HERNQUIST', 'PJAFFE', 'PJAFFE_ELLIPSE',
+            if model in ['SERSIC', 'SERSIC_ELLIPSE', 'SERSIC_ELLIPSE_Q_PHI', 'CORE_SERSIC',
+                         'HERNQUIST', 'PJAFFE', 'PJAFFE_ELLIPSE',
                          'HERNQUIST_ELLIPSE', 'GAUSSIAN', 'GAUSSIAN_ELLIPSE', 'POWER_LAW', 'NIE', 'CHAMELEON',
                          'DOUBLE_CHAMELEON', 'TRIPLE_CHAMELEON', 'UNIFORM', 'INTERPOL', 'ELLIPSOID']:
                 kwargs_list[k]['amp'] = param[i]
diff --git a/lenstronomy/Sampling/Likelihoods/SKiNN_1.0_config.json b/lenstronomy/Sampling/Likelihoods/SKiNN_1.0_config.json
new file mode 100644
index 000000000..e8617b34e
--- /dev/null
+++ b/lenstronomy/Sampling/Likelihoods/SKiNN_1.0_config.json
@@ -0,0 +1,1 @@
+{"grid_settings": {"deltapix": 0.02, "npix": 551}, "training_abs_bounds": {"q_mass": [0.6, 1.0], "q_light": [0.6, 1.0], "theta_E": [0.5, 2.0], "n_sersic": [2.0, 4.0], "R_sersic": [0.25, 2], "core_size": [0.0, 0.081], "gamma": [0.25, 0.75], "b_ani": [-0.4, 0.4], "incli": [0, 90]}, "other_notes": {"PAs": "PA of mass and light must be aligned", "R_sersic": "Must be between 0.5*theta_E and theta_E", "incli": "Must be greater than arccos(q) for q_light and q_mass"}}
\ No newline at end of file
diff --git a/lenstronomy/Sampling/Likelihoods/kinematic_2D_likelihood.py b/lenstronomy/Sampling/Likelihoods/kinematic_2D_likelihood.py
new file mode 100644
index 000000000..be585f66c
--- /dev/null
+++ b/lenstronomy/Sampling/Likelihoods/kinematic_2D_likelihood.py
@@ -0,0 +1,230 @@
+import numpy as np
+import lenstronomy.Util.param_util as param_util
+import lenstronomy.Util.util as util
+from scipy import signal
+from lenstronomy.Util.kin_sampling_util import KinNNImageAlign
+from lenstronomy.Sampling.Likelihoods import kinematic_NN_call
+
+__all__ = ['KinLikelihood']
+
+
+class KinLikelihood(object):
+    """
+    Class to compute the likelihood associated to binned 2D kinematic maps
+    """
+
+    def __init__(self, kinematic_bin_2d_class, lens_model_class, lens_light_model_class, kwargs_data, idx_lens=0,
+                 idx_lens_light=0):
+        """
+        :param kinematic_data_2d_class: KinBin class instance
+        :param lens_model_class: LensModel class instance
+        :param lens_light_model_class: LightModel class instance
+        :param kwargs_data: kwargs describing image rotation
+        :param idx_lens: int, index of the LensModel mass profile to consider for kinematics
+        :param idx_lens_light: int, index of the lens LightModel profile to consider for kinematics
+        """
+        self.lens_model_class = lens_model_class
+        self.z_lens = self.lens_model_class.z_lens
+        self.lens_light_model_class = lens_light_model_class
+        self._idx_lens = idx_lens
+        self._idx_lens_light = idx_lens_light
+        self.kin_class = kinematic_bin_2d_class
+
+        # setup KinematicNN and KinNNImagealign class inputs to be updated when called
+        self.kin_input = self.kin_class.kin_bin2kwargs()
+        self.image_input = self.kwargs_data2image_input(kwargs_data)
+        self.kinematic_NN = kinematic_NN_call.KinematicNN()
+        self.config = self.kinematic_NN.config
+        self.kinNN_input = {'deltaPix': self.config['grid_settings']['deltapix'],
+                            'image': np.ones((self.config['grid_settings']['npix'], self.config['grid_settings']['npix']))}
+        self.KiNNalign = KinNNImageAlign(self.kin_input, self.image_input, self.kinNN_input)
+
+        self.data = self.kin_class.data
+        self.psf = self.kin_class.PSF.kernel_point_source
+        self.bin_mask = self.kin_class.bin_mask
+        self.covariance = self.kin_class.covariance
+
+        self.kin_x_grid, self.kin_y_grid = self.kin_class.kin_grid()
+
+        self.lens_light_bool_list = list(np.zeros_like(self.lens_light_model_class.profile_type_list, dtype=bool))
+        self.lens_light_bool_list[self._idx_lens_light] = True
+
+        d_d_fiducial = 1215.739  # fiducial distances used to create training set
+        d_s_fiducial = 1650.753
+        d_ds_fiducial = 1042.883
+        z_d_fiducial = 0.5
+        d_dt_fiducial = (1 + z_d_fiducial) * d_d_fiducial * d_s_fiducial / d_ds_fiducial
+        self.fiducial_scale = d_dt_fiducial / (d_d_fiducial * (1 + z_d_fiducial))
+        self.vrms = None
+
+    def calc_vrms(self, kwargs_lens, kwargs_lens_light, kwargs_special, verbose=False):
+        """
+        Calculates binned vrms using SKiNN
+
+        :param kwargs_lens: lens model kwargs list
+        :param kwargs_lens_light: lens light kwargs list
+        :param kwargs_special: cosmology and other kwargs
+        :param verbose: default False; if True print statements when out of bounds
+        return binned vrms [km/s]; if SKiNN not installed return nan
+        """
+        self.update_image_input(kwargs_lens)
+        self.light_map = self.lens_light_model_class.surface_brightness(self.kin_x_grid, self.kin_y_grid,
+                                                                        kwargs_lens_light,
+                                                                        self.lens_light_bool_list)
+        input_params, same_orientation = self.convert_to_nn_params(kwargs_lens, kwargs_lens_light, kwargs_special)
+        if self.kinematic_NN.SKiNN_installed:
+            velo_map = self.kinematic_NN.generate_map(input_params, verbose=verbose)
+            velo_map = self.rescale_distance(velo_map, kwargs_special)  # RESCALE ACCORDING TO D_d, D_dt
+            # Rotation and interpolation in kin data coordinates
+            self.kinNN_input['image'] = velo_map
+            self.KiNNalign.update(self.kin_input, self.image_input, self.kinNN_input)
+            self.rotated_velo = self.KiNNalign.interp_image()
+            # Convolution by PSF to calculate Vrms and binning
+            vrms = self.auto_binning(self.rotated_velo, self.light_map)
+            return vrms
+        else:
+            return np.nan
+
+    def logL(self, kwargs_lens, kwargs_lens_light, kwargs_special, verbose=False):
+        """
+        Calculates Log likelihood from 2D kinematic likelihood
+
+        :param kwargs_lens: lens model kwargs list
+        :param kwargs_lens_light: lens light kwargs list
+        :param kwargs_special: cosmology and other kwargs
+        :param verbose: default False; if True print statements when out of bounds
+        return kinematics log likelihood
+        """
+        if self.kinematic_NN.SKiNN_installed:
+            self.update_image_input(kwargs_lens)
+            self.light_map = self.lens_light_model_class.surface_brightness(self.kin_x_grid, self.kin_y_grid,
+                                                                            kwargs_lens_light,
+                                                                            self.lens_light_bool_list)
+            input_params, same_orientation = self.convert_to_nn_params(kwargs_lens, kwargs_lens_light, kwargs_special)
+            if not self.kinematic_NN.check_bounds(input_params, same_orientation=same_orientation,
+                                                  verbose=verbose):
+                # params not within training set. Penalty
+                return -10 ** 8
+            self.vrms = self.calc_vrms(kwargs_lens, kwargs_lens_light, kwargs_special, verbose=verbose)
+            logL = self._logL(self.vrms)
+        else:
+            logL = np.nan
+        return logL
+
+    def convert_to_nn_params(self, kwargs_lens, kwargs_lens_light, kwargs_special):
+        """
+        converts lenstronomy kwargs into input vector for SKiNN, also returns whether or not mass and light are aligned
+
+        :param kwargs_lens: lens model kwargs list
+        :param kwargs_lens_light: lens light kwargs list
+        :param kwargs_special: cosmology and other kwargs
+        :return: parameters in GLEE convention to be input into NN
+        """
+        # lenstronomy to GLEE conversion
+        # orientation_mass,q_mass=param_util.ellipticity2phi_q(kwargs_lens[self._idx_lens]['e1'],
+        #                                                      kwargs_lens[self._idx_lens]['e2'])
+        orientation_mass = kwargs_lens[self._idx_lens]['phi']
+        q_mass = kwargs_lens[self._idx_lens]['q']
+        # orientation_light, q_light = param_util.ellipticity2phi_q(kwargs_lens_light[self._idx_lens]['e1'],
+        #                                                         kwargs_lens_light[self._idx_lens]['e2'])
+        orientation_light = kwargs_lens_light[self._idx_lens]['phi']
+        q_light = kwargs_lens_light[self._idx_lens]['q']
+
+        same_orientation = True  # confirm that angles are aligned
+        if np.abs(orientation_mass - orientation_light) > 0.05:
+            if q_light < 0.95 or q_mass < 0.95:
+                same_orientation = False
+
+        theta_e_lenstro = kwargs_lens[self._idx_lens]['theta_E']
+        gamma_lenstro = kwargs_lens[self._idx_lens]['gamma']
+        gamma_glee = (gamma_lenstro - 1) / 2
+        r_e_scale = (2 / (1 + q_mass)) ** (1 / (2 * gamma_glee)) * np.sqrt(q_mass)
+        theta_e_glee = theta_e_lenstro / r_e_scale
+        return np.array([q_mass, q_light, theta_e_glee, kwargs_lens_light[self._idx_lens]['n_sersic'],
+                         kwargs_lens_light[self._idx_lens]['R_sersic'], 8.0e-2, gamma_glee,
+                         kwargs_special['b_ani'], kwargs_special['incli'] * 180 / np.pi]), same_orientation
+
+    def rescale_distance(self, image, kwargs_special):
+        """
+        rescales velocity map according to distance, requires lens redshift
+
+        :param image: vrms image [km/s]
+        :param kwargs_special: kwargs with cosmological distances for rescaling
+        return rescaled vrms image [km/s]
+        """
+        new_scale = kwargs_special['D_dt'] / (kwargs_special['D_d'] * (1 + self.z_lens))
+        factor = np.sqrt(
+            new_scale / self.fiducial_scale)  # vrms^2 propto D_dt/Dd (see https://arxiv.org/pdf/2109.14615.pdf Eq. 20)
+        return image * factor
+
+    def kwargs_data2image_input(self, kwargs_data):
+        """
+        Creates the kwargs of the image needed for 2D kinematic likelihood
+
+        :param kwargs_data: kwargs giving image and describing imaging data coordinate transformation
+        :return kwargs: coordinate transformation kwargs as input for KinNNImageAlign class
+        """
+        delta_pix = np.sqrt(np.abs(np.linalg.det(kwargs_data['transform_pix2angle'])))
+        kwargs = {'image': kwargs_data['image_data'], 'deltaPix': delta_pix,
+                  'transform_pix2angle': kwargs_data['transform_pix2angle'],
+                  'ra_at_xy0': kwargs_data['ra_at_xy_0'], 'dec_at_xy0': kwargs_data['dec_at_xy_0']}
+        return kwargs
+
+    def update_image_input(self, kwargs_lens):
+        """
+        Updates the image_input for rotation with the new values of orientation and center of the lens model.
+
+        :param kwargs_lens: lens kwargs with center and PA positions to be used as input
+        """
+        orientation_ellipse = kwargs_lens[self._idx_lens]['phi']
+        q = kwargs_lens[self._idx_lens]['q']
+        cx = kwargs_lens[self._idx_lens]['center_x']
+        cy = kwargs_lens[self._idx_lens]['center_y']
+        self.image_input['ellipse_PA'] = orientation_ellipse
+        self.image_input['offset_x'] = cx
+        self.image_input['offset_y'] = cy
+
+    def auto_binning(self, rotated_map, light_map):
+        """
+        Function to convolve and bin the NN rotated output
+
+        :param rotated_map: model vrms map [km/s] in data pixel coordinates
+        :param light_map: model light map in data pixel coordinates for weighting
+        :return: binned vrms [km/s] for comparison with data
+        """
+        vrms = rotated_map
+        mge_car = light_map
+
+        vrms = signal.fftconvolve(vrms, self.psf, mode='same')
+        mge_car_con = signal.fftconvolve(mge_car, self.psf, mode='same')
+
+        numerator = []
+        denominator = []
+
+        for idx in range(len(self.data)):
+            math_pos = np.where(self.bin_mask == idx)
+            if np.shape(math_pos[0])[0] == 0:
+                raise ValueError('binmap mismatch with data: no pixels in bin with idx %i' % (idx))
+            numerator.append(np.sum(vrms[math_pos] * mge_car_con[math_pos]))
+            denominator.append(np.sum(mge_car_con[math_pos]))
+
+        vrms = np.array(numerator) / np.array(denominator).clip(0)
+
+        return vrms
+
+    def _logL(self, vrms):
+        """
+        Calculates the log likelihood for a given binned model
+
+        :param vrms: binned vrms [km/s] to compare with observed binned data
+        :return: log likelihood
+        """
+        # log_like = (vrms - self.data)**2 / self.noise**2
+        cov_inv = np.linalg.inv(self.covariance)
+        log_like = np.matmul(np.matmul((vrms - self.data).T, cov_inv), (vrms - self.data))
+        logL = - np.sum(log_like) / 2
+
+        if not np.isfinite(logL):
+            return -10 ** 15
+
+        return logL
diff --git a/lenstronomy/Sampling/Likelihoods/kinematic_NN_call.py b/lenstronomy/Sampling/Likelihoods/kinematic_NN_call.py
new file mode 100644
index 000000000..f42b6cfbb
--- /dev/null
+++ b/lenstronomy/Sampling/Likelihoods/kinematic_NN_call.py
@@ -0,0 +1,77 @@
+# import trained NN
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+import os
+import json
+
+
+class KinematicNN():
+    """
+    Class to call the NN to emulate JAM kinematics
+    """
+
+    def __init__(self):
+        current_directory = os.path.dirname(os.path.abspath(__file__))
+        filepath = os.path.join(current_directory, "SKiNN_1.0_config.json")
+        with open(filepath, 'r') as json_file:
+            self.config = json.load(json_file)
+        try:
+            import SKiNN
+            from SKiNN.generator import Generator
+            self.SKiNN_installed = True
+            self.generator = Generator()
+        except:
+            print("Warning : SKiNN not installed properly, \
+        but tests will be trivially fulfilled. \
+                  Get SKiNN from https://github.com/mattgomer/SKiNN")
+            self.SKiNN_installed = False
+
+    def generate_map(self, input_p, verbose=False):
+        """Generate velocity map given input parameters.
+
+        :param input_p: vector of input parameters
+        :return: the velocity map associated with input_p
+        """
+        self.within_bounds = self.check_bounds(input_p, verbose)
+        return self.generator.generate_map(input_p)
+
+    def plot_map(self, ax, input_p):
+        self.check_bounds(input_p)
+        im = ax.imshow(self.generate_map(input_p))
+        ax.set_title('Prediction')
+        divider = make_axes_locatable(ax)
+        cax = divider.append_axes("right", size="5%", pad=0.05)
+        cb = plt.colorbar(im, cax=cax)
+        cb.set_label(r'$v_{\rm rms}$ [km/s]')
+        return ax
+
+    def check_bounds(self, input_p, same_orientation=True, verbose=False):
+        """
+        Checks to see if input parameters lie in bounds used for the training set
+
+        :param input_p: input parameters to NN
+        :param same_orientation: default True; confirms that mass and light have same position angles
+        :param verbose: default False; if True prints statements when out of bounds
+        """
+        within_bounds = True
+        if not same_orientation:
+            if verbose:
+                print('NN CALL WARNING: Mass and light have different PAs!')
+            within_bounds = False
+        training_abs_bounds=self.config['training_abs_bounds']
+        for idx, pname in enumerate(['q_mass', 'q_light', 'theta_E', 'n_sersic', 'R_sersic',
+                                     'core_size', 'gamma', 'b_ani', 'incli']):
+            if input_p[idx] < training_abs_bounds[pname][0] or input_p[idx] > training_abs_bounds[pname][1]:
+                if verbose:
+                    print('NN CALL WARNING: param', pname, 'is outside of training bounds!')
+                within_bounds = False
+        if input_p[4] > input_p[2] or input_p[4] < 0.5 * input_p[2]:
+            if verbose:
+                print('NN CALL WARNING: R_sersic is not within training bounds- must be between 0.5R_E and R_E!')
+            within_bounds = False
+        if (input_p[8] - 2) * np.pi / 180 < np.arccos(np.min([input_p[0], input_p[1]])):  # wrt q_light and q_mass
+            if verbose:
+                print('NN CALL WARNING: Inclination is nonphysical!')
+            within_bounds = False
+        return within_bounds
diff --git a/lenstronomy/Sampling/likelihood.py b/lenstronomy/Sampling/likelihood.py
index 575c12e70..670aa152f 100644
--- a/lenstronomy/Sampling/likelihood.py
+++ b/lenstronomy/Sampling/likelihood.py
@@ -5,6 +5,7 @@
 from lenstronomy.Sampling.Likelihoods.position_likelihood import PositionLikelihood
 from lenstronomy.Sampling.Likelihoods.flux_ratio_likelihood import FluxRatioLikelihood
 from lenstronomy.Sampling.Likelihoods.prior_likelihood import PriorLikelihood
+from lenstronomy.Sampling.Likelihoods.kinematic_2D_likelihood import KinLikelihood
 import lenstronomy.Util.class_creator as class_creator
 import numpy as np
 
@@ -22,6 +23,7 @@ class LikelihoodModule(object):
 
     Additional arguments are supported for adding a time-delay likelihood etc (see __init__ definition)
     """
+
     def __init__(self, kwargs_data_joint, kwargs_model, param_class, image_likelihood=True, check_bounds=True,
                  check_matched_source_position=False, astrometric_likelihood=False, image_position_likelihood=False,
                  source_position_likelihood=False, image_position_uncertainty=0.004, check_positive_flux=False,
@@ -34,7 +36,8 @@ def __init__(self, kwargs_data_joint, kwargs_model, param_class, image_likelihoo
                  prior_source_kde=None, prior_lens_light_kde=None, prior_ps_kde=None, prior_special_kde=None,
                  prior_extinction_kde=None, prior_lens_lognormal=None, prior_source_lognormal=None,
                  prior_extinction_lognormal=None, prior_lens_light_lognormal=None, prior_ps_lognormal=None,
-                 prior_special_lognormal=None, custom_logL_addition=None, kwargs_pixelbased=None):
+                 prior_special_lognormal=None, custom_logL_addition=None, kwargs_pixelbased=None,
+                 kinematic_2d_likelihood=False, kin_lens_idx=0, kin_lens_light_idx=0):
         """
         initializing class
 
@@ -75,10 +78,13 @@ def __init__(self, kwargs_data_joint, kwargs_model, param_class, image_likelihoo
          kwargs_ps, kwargs_special, kwargs_extinction) and returns a logL (punishing) value.
         :param kwargs_pixelbased: keyword arguments with various settings related to the pixel-based solver
          (see SLITronomy documentation)
+        :param kinematic_2d_likelihood: bool, option to compute the kinematic likelihood
         """
-        multi_band_list, multi_band_type, time_delays_measured, time_delays_uncertainties, flux_ratios, flux_ratio_errors, ra_image_list, dec_image_list = self._unpack_data(**kwargs_data_joint)
+        multi_band_list, multi_band_type, time_delays_measured, time_delays_uncertainties, flux_ratios, flux_ratio_errors, ra_image_list, dec_image_list, kinematic_data = self._unpack_data(
+            **kwargs_data_joint)
         if len(multi_band_list) == 0:
             image_likelihood = False
+        self.kinematic_data = kinematic_data
         self.param = param_class
         self._lower_limit, self._upper_limit = self.param.param_limits()
         self._prior_likelihood = PriorLikelihood(prior_lens, prior_source, prior_lens_light, prior_ps, prior_special,
@@ -92,6 +98,7 @@ def __init__(self, kwargs_data_joint, kwargs_model, param_class, image_likelihoo
         self._time_delay_likelihood = time_delay_likelihood
         self._image_likelihood = image_likelihood
         self._flux_ratio_likelihood = flux_ratio_likelihood
+        self._kinematic_2D_likelihood = kinematic_2d_likelihood
         if kwargs_flux_compute is None:
             kwargs_flux_compute = {}
         linear_solver = self.param.linear_solver
@@ -117,11 +124,25 @@ def __init__(self, kwargs_data_joint, kwargs_model, param_class, image_likelihoo
                                  'restrict_image_number': restrict_image_number, 'max_num_images': max_num_images}
         self._kwargs_flux = {'flux_ratios': flux_ratios, 'flux_ratio_errors': flux_ratio_errors}
         self._kwargs_flux.update(self._kwargs_flux_compute)
+
+        if self._kinematic_2D_likelihood is True:
+            print(
+                "Note that the 2D kinematic likelihood assumes that the lens and lens light have the same center and orientation")
+            if len(multi_band_list) > 1:
+                print('Kinematic Likelihood not meant for multiband, using first band by default')
+            if kwargs_model['lens_model_list'][kin_lens_idx] not in ['EPL_Q_PHI']:
+                print('Lens for kinematic is not EPL_Q_PHI, the 2D kinematic likelihood will break.')
+            if kwargs_model['lens_light_model_list'][kin_lens_light_idx] not in ['SERSIC_ELLIPSE_Q_PHI']:
+                print('Lens light for kinematic is not SERSIC_ELLIPSE_Q_PHI , the 2D kinematic likelihood will break.')
+            self._kin_lens_idx = kin_lens_idx
+            self._kin_lens_light_idx = kin_lens_light_idx
+
         self._class_instances(kwargs_model=kwargs_model, kwargs_imaging=self.kwargs_imaging,
                               kwargs_position=self._kwargs_position, kwargs_flux=self._kwargs_flux,
-                              kwargs_time_delay=self._kwargs_time_delay)
+                              kwargs_time_delay=self._kwargs_time_delay, kinematic_data=self.kinematic_data)
 
-    def _class_instances(self, kwargs_model, kwargs_imaging, kwargs_position, kwargs_flux, kwargs_time_delay):
+    def _class_instances(self, kwargs_model, kwargs_imaging, kwargs_position, kwargs_flux, kwargs_time_delay,
+                         kinematic_data):
         """
 
         :param kwargs_model: lenstronomy model keyword arguments
@@ -129,13 +150,15 @@ def _class_instances(self, kwargs_model, kwargs_imaging, kwargs_position, kwargs
         :param kwargs_position: keyword arguments for positional likelihood
         :param kwargs_flux: keyword arguments for flux ratio likelihood
         :param kwargs_time_delay: keyword arguments for time delay likelihood
+        :param kinematic_data: kinematic class for kinematic likelihood
         :return: updated model instances of this class
         """
 
         # TODO: in case lens model or point source models are only applied on partial images, then this current class
         # has ambiguities when it comes to position likelihood, time-delay likelihood and flux ratio likelihood
-        lens_model_class, _, _, point_source_class, _ = class_creator.create_class_instances(all_models=True,
-                                                                                             **kwargs_model)
+        lens_model_class, _, lens_light_model_class, point_source_class, _ = class_creator.create_class_instances(
+            all_models=True,
+            **kwargs_model)
         self.PointSource = point_source_class
 
         if self._time_delay_likelihood is True:
@@ -148,6 +171,10 @@ def _class_instances(self, kwargs_model, kwargs_imaging, kwargs_position, kwargs
         self._position_likelihood = PositionLikelihood(point_source_class, **kwargs_position)
         if self._flux_ratio_likelihood is True:
             self.flux_ratio_likelihood = FluxRatioLikelihood(lens_model_class, **kwargs_flux)
+        if self._kinematic_2D_likelihood is True:
+            self.kinematic_2D_likelihood = KinLikelihood(kinematic_data, lens_model_class, lens_light_model_class,
+                                                         kwargs_imaging['multi_band_list'][0][0], self._kin_lens_idx,
+                                                         self._kin_lens_light_idx)
 
     def __call__(self, a):
         return self.logL(a)
@@ -168,7 +195,7 @@ def logL(self, args, verbose=False):
         if self._check_bounds is True:
             penalty, bound_hit = self.check_bounds(args, self._lower_limit, self._upper_limit, verbose=verbose)
             if bound_hit is True:
-                return -10**15
+                return -10 ** 15
         return self.log_likelihood(kwargs_return, verbose=verbose)
 
     def log_likelihood(self, kwargs_return, verbose=False):
@@ -186,10 +213,10 @@ def log_likelihood(self, kwargs_return, verbose=False):
          - 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'], \
-                                                                                   kwargs_return['kwargs_ps'], \
-                                                                                   kwargs_return['kwargs_special']
+            kwargs_return['kwargs_source'], \
+            kwargs_return['kwargs_lens_light'], \
+            kwargs_return['kwargs_ps'], \
+            kwargs_return['kwargs_special']
         # update model instance in case of changes affecting it (i.e. redshift sampling in multi-plane)
         self._update_model(kwargs_special)
         # generate image and computes likelihood
@@ -214,6 +241,11 @@ def log_likelihood(self, kwargs_return, verbose=False):
             logL += logL_flux_ratios
             if verbose is True:
                 print('flux ratio logL = %s' % logL_flux_ratios)
+        if self._kinematic_2D_likelihood is True:
+            logL_kinematic_2d = self.kinematic_2D_likelihood.logL(kwargs_lens, kwargs_lens_light, kwargs_special)
+            logL += logL_kinematic_2d
+            if verbose is True:
+                print('kinematic logL = %s' % logL_kinematic_2d)
         logL += self._position_likelihood.logL(kwargs_lens, kwargs_ps, kwargs_special, verbose=verbose)
         logL_prior = self._prior_likelihood.logL(**kwargs_return)
         logL += logL_prior
@@ -237,7 +269,7 @@ def check_bounds(args, lowerLimit, upperLimit, verbose=False):
         args = np.atleast_1d(args)
         for i in range(0, len(args)):
             if args[i] < lowerLimit[i] or args[i] > upperLimit[i]:
-                penalty = 10.**5
+                penalty = 10. ** 5
                 bound_hit = True
                 if verbose is True:
                     print('parameter %s with value %s hit the bounds [%s, %s] ' % (i, args[i], lowerLimit[i],
@@ -290,7 +322,7 @@ def negativelogL(self, a):
     @staticmethod
     def _unpack_data(multi_band_list=None, multi_band_type='multi-linear', time_delays_measured=None,
                      time_delays_uncertainties=None, flux_ratios=None, flux_ratio_errors=None, ra_image_list=None,
-                     dec_image_list=None):
+                     dec_image_list=None, kinematic_data=None):
         """
 
         :param multi_band_list: list of [[kwargs_data, kwargs_psf, kwargs_numerics], [], ...]
@@ -308,7 +340,7 @@ def _unpack_data(multi_band_list=None, multi_band_type='multi-linear', time_dela
         if dec_image_list is None:
             dec_image_list = []
         return multi_band_list, multi_band_type, time_delays_measured, time_delays_uncertainties, flux_ratios, \
-               flux_ratio_errors, ra_image_list, dec_image_list
+            flux_ratio_errors, ra_image_list, dec_image_list, kinematic_data
 
     def _reset_point_source_cache(self, bool_input=True):
         self.PointSource.delete_lens_model_cache()
@@ -322,11 +354,11 @@ def _update_model(self, kwargs_special):
         modeled redshifts of the deflector and/or source planes are made
 
         :param kwargs_special: keyword arguments from SpecialParam() class return of sampling arguments
-        :return: None, all class instances updated to recent modek
+        :return: None, all class instances updated to recent model
         """
         kwargs_model, update_bool = self.param.update_kwargs_model(kwargs_special)
         if update_bool is True:
             self._class_instances(kwargs_model=kwargs_model, kwargs_imaging=self.kwargs_imaging,
                                   kwargs_position=self._kwargs_position, kwargs_flux=self._kwargs_flux,
-                                  kwargs_time_delay=self._kwargs_time_delay)
+                                  kwargs_time_delay=self._kwargs_time_delay, kinematic_data=self.kinematic_data)
         # TODO remove redundancies with Param() calls updates
diff --git a/lenstronomy/Sampling/parameters.py b/lenstronomy/Sampling/parameters.py
index 398b15b33..9dbb7c239 100644
--- a/lenstronomy/Sampling/parameters.py
+++ b/lenstronomy/Sampling/parameters.py
@@ -122,7 +122,7 @@ def __init__(self, kwargs_model,
                  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,
+                 source_redshift_sampling_indexes=None, source_grid_offset=False, kinematic_sampling=None, num_shapelet_lens=0,
                  log_sampling_lens=[]):
         """
 
@@ -331,7 +331,7 @@ def __init__(self, kwargs_model,
                                           kwargs_lower=kwargs_lower_special, kwargs_upper=kwargs_upper_special,
                                           point_source_offset=self._point_source_offset, num_images=self._num_images,
                                           source_size=source_size, num_tau0=num_tau0, num_z_sampling=num_z_sampling,
-                                          source_grid_offset=source_grid_offset)
+                                          source_grid_offset=source_grid_offset,kinematic_sampling=kinematic_sampling)
         for lens_source_joint in self._joint_lens_with_source_light:
             i_source = lens_source_joint[0]
             if i_source in self._image_plane_source_list:
diff --git a/lenstronomy/Sampling/special_param.py b/lenstronomy/Sampling/special_param.py
index c8841e4ea..a4c7c84a2 100644
--- a/lenstronomy/Sampling/special_param.py
+++ b/lenstronomy/Sampling/special_param.py
@@ -19,6 +19,29 @@ class DdtSamplingParam(SingleParam):
     _kwargs_lower = {'D_dt': 0}
     _kwargs_upper = {'D_dt': 100000}
 
+class DdSamplingParam(SingleParam):
+    '''
+    Deflector distance parameter
+    '''
+    param_names = ['D_d']
+    _kwargs_lower = {'D_d': 0}
+    _kwargs_upper = {'D_d': 100000}
+
+class BetaAnisotropyParam(SingleParam):
+    '''
+    Cylindrical anisotropy parameter
+    '''
+    param_names = ['b_ani']
+    _kwargs_lower = {'b_ani': -1}
+    _kwargs_upper = {'b_ani': 1}
+
+class InclinationParam(SingleParam):
+    '''
+    Inclination parameter (radians)
+    '''
+    param_names = ['incli']
+    _kwargs_lower = {'incli': 0}
+    _kwargs_upper = {'incli': np.pi/2}
 
 class SourceSizeParam(SingleParam):
     '''
@@ -139,7 +162,7 @@ class that handles special parameters that are not directly part of a specific m
     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):
+                 num_z_sampling=0, source_grid_offset=False, kinematic_sampling=False):
         """
 
         :param Ddt_sampling: bool, if True, samples the time-delay distance D_dt (in units of Mpc)
@@ -158,9 +181,16 @@ def __init__(self, Ddt_sampling=False, mass_scaling=False, num_scale_factor=1,
         :param source_grid_offset: bool, if True, samples two parameters (x, y) for the offset of the pixelated source
          plane grid coordinates.
          Warning: this is only defined for pixel-based source modelling (e.g. 'SLIT_STARLETS' light profile)
+        :param kinematic_sampling: bool, if True, samples the kinematic parameters b_ani, incli, with cosmography
+         D_dt (overrides _D_dt_sampling) and Dd
         """
 
-        self._D_dt_sampling = DdtSamplingParam(Ddt_sampling)
+
+        self._D_dt_sampling = DdtSamplingParam(Ddt_sampling or kinematic_sampling)
+
+        self._D_d_sampling = DdSamplingParam(kinematic_sampling)
+        self._b_ani_sampling = BetaAnisotropyParam(kinematic_sampling)
+        self._incli_sampling = InclinationParam(kinematic_sampling)
         if not mass_scaling:
             num_scale_factor = 0
         self._mass_scaling = MassScalingParam(num_scale_factor)
@@ -233,6 +263,9 @@ def num_param(self):
     @property
     def _param_groups(self):
         return [self._D_dt_sampling,
+                self._D_d_sampling,
+                self._b_ani_sampling,
+                self._incli_sampling,
                 self._mass_scaling,
                 self._general_scaling,
                 self._point_source_offset,
diff --git a/lenstronomy/Util/image_util.py b/lenstronomy/Util/image_util.py
index 77d6d4cfd..f449b3c1b 100644
--- a/lenstronomy/Util/image_util.py
+++ b/lenstronomy/Util/image_util.py
@@ -117,6 +117,17 @@ def rotateImage(img, angle):
     imgR = ndimage.rotate(img, angle, reshape=False)
     return imgR
 
+@export
+def shift_image(img, shift):
+    """
+    queries scipy.ndimage.shift routine
+
+    :param img: image to be shifted
+    :param shift: sequence containing x and y shift in pixels
+    :return: shifted image
+    """
+    img_s = ndimage.shift(img, shift)
+    return img_s
 
 @export
 def re_size_array(x_in, y_in, input_values, x_out, y_out):
diff --git a/lenstronomy/Util/kin_sampling_util.py b/lenstronomy/Util/kin_sampling_util.py
new file mode 100644
index 000000000..219cb7dfa
--- /dev/null
+++ b/lenstronomy/Util/kin_sampling_util.py
@@ -0,0 +1,220 @@
+__author__ = 'Matt Gomer'
+
+import numpy as np
+from scipy.interpolate import RectBivariateSpline
+import matplotlib.pyplot as plt
+
+
+class KinNNImageAlign(object):
+    """
+    class to rotate and interpolate SKiNN image (which is built along x-axis) to be
+        on the grid of the spectral data (e.g., MUSE), rotated to match light distribution as defined by the lens model
+    main function is interp_image() which will output a 2D image interpolated on the spectra grid
+    """
+
+    def __init__(self, spectra_inputs, imaging_inputs, kin_nn_inputs):
+        """
+        initialize input data
+
+        :param spectra_inputs: dictionary which encodes grid and transformation information for kinematic data
+            :'image': contains 2d image used to calculate grid coordinates
+            :'transform_pix2angle': transformation matrix to convert from pixel xy to ra/dec
+            :'ra_at_xy0': ra coordinate at pixel (0,0)
+            :'dec_at_xy0': dec coordinate at pixel (0,0)
+
+        :param imaging_inputs: dictionary which encodes grid and transformation information for imaging data
+            :'image': contains 2d image used to calculate grid coordinates
+            :'transform_pix2angle': transformation matrix to convert from pixel xy to ra/dec
+            :'ra_at_xy0': ra coordinate at pixel (0,0)
+            :'dec_at_xy0': dec coordinate at pixel (0,0)
+            :'ellipse_PA': position angle of ellipse axis relative to x direction
+            :'offset_x': how many pixels to offset the center of the grid to match the kinNN center (x-direction)
+            :'offset_y': how many pixels to offset the center of the grid to match the kinNN center (y-direction)
+
+        :param kin_nn_inputs: dictionary which encodes grid information for NN output data
+            :'image': contains 2d image used to calculate grid coordinates
+            :'deltaPix': pixel size
+
+        """
+        self.spectra_data = spectra_inputs
+        self.imaging_data = imaging_inputs
+        self.imaging_deltapix = np.sqrt(np.abs(np.linalg.det(self.imaging_data['transform_pix2angle'])))
+        self.kinNN_data = kin_nn_inputs
+        self.write_npix()
+
+    def update(self, spectra_inputs=None, imaging_inputs=None, kin_nn_inputs=None, update_npix=False):
+        """
+        Update with inputs
+        """
+        if spectra_inputs is not None:
+            self.spectra_data = spectra_inputs
+        if imaging_inputs is not None:
+            self.imaging_data = imaging_inputs
+        if kin_nn_inputs is not None:
+            self.kinNN_data = kin_nn_inputs
+        if update_npix:
+            self.write_npix()
+
+    def write_npix(self):
+        """
+        Check that images are squared and write the keyword 'npix'
+        """
+        for input_set in [self.spectra_data, self.imaging_data, self.kinNN_data]:
+            # make sure each image is square and add npix to each dictionary
+            if 'image' in input_set.keys():
+                if np.shape(input_set['image'])[0] != np.shape(input_set['image'])[1]:
+                    raise ValueError('current version only works for square images')
+                npix = np.shape(input_set['image'])[0]
+                input_set['npix'] = npix
+
+    def pix_coords(self, input_set, flatten=True):
+        """
+        simple function to give pixel coordinates of grid
+
+        :param input_set: dictionary from above (e.g. spectra_inputs) which completely describes grid transformation
+        :boolean flatten: default True; if True, return 1D flattened output, if False, return 2D grid
+        :return pixel coordinates of grid
+        """
+        x_grid = np.tile(np.arange(input_set['npix']), input_set['npix'])
+        y_grid = np.repeat(np.arange(input_set['npix']), input_set['npix'])
+        if flatten is True:
+            return x_grid, y_grid
+        else:
+            return x_grid.reshape(input_set['npix'], input_set['npix']), y_grid.reshape(input_set['npix'],
+                                                                                        input_set['npix'])
+
+    def radec_to_xy(self, ra, dec, xy_to_radec_matrix, ra_atxy0, dec_atxy0):
+        """
+        converts from radec to pixel coordinates
+
+        :param ra: ra coordinate to transform
+        :param dec: dec coordinate to transform
+        :param xy_to_radec_matrix: transformation matrix to convert from pixel xy to ra/dec
+        :param ra_atxy0: ra coordinate at pixel (0,0)
+        :param dec_atxy0: dec coordinate at pixel (0,0)
+
+        :return: x and y coordinates
+        """
+        ra = ra - ra_atxy0
+        dec = dec - dec_atxy0
+        x, y = np.linalg.inv(xy_to_radec_matrix).dot(np.array([ra, dec]))
+        return x, y
+
+    def xy_to_radec(self, x, y, xy_to_radec_matrix, ra_atxy0, dec_atxy0):
+        """
+        converts from pixel coordinates to radec
+
+        :param x: x coordinate to transform
+        :param y: y coordinate to transform
+        :param xy_to_radec_matrix: transformation matrix to convert from pixel xy to ra/dec
+        :param ra_atxy0: ra coordinate at pixel (0,0)
+        :param dec_atxy0: dec coordinate at pixel (0,0)
+
+        :return: ra and dec coordinates
+        """
+        ra, dec = xy_to_radec_matrix.dot(np.array([x, y]))
+        return ra + ra_atxy0, dec + dec_atxy0
+
+    def rotate_imaging_into_kin_nn(self, imaging_x, imaging_y, ellipse_pa_to_imagingx_angle,
+                                   deltapix_imaging, deltapix_kin_nn, npix_imaging, npix_kin_nn,
+                                   offsetx=0, offsety=0):
+        """
+        rotates and rescales from the x,y imaging coordinate system into the NN coordinate system
+
+        :param imaging_x: imaging x coordinate to transform
+        :param imaging_y: imaging y coordinate to transform
+        :param ellipse_pa_to_imagingx_angle: (radians) position angle of ellipse major axis relative to x
+            in the imaging coordinate system
+        :param deltapix_imaging: pixel size of imaging image
+        :param deltapix_kin_nn: pixel size of NN image
+        :param npix_imaging: number of pixels on a side of the imaging image
+        :param npix_kin_nn: number of pixels on a side of the NN image
+        :param offsetx: how many pixels to offset the center of the grid to match the kinNN center (x-direction)
+        :param offsety: how many pixels to offset the center of the grid to match the kinNN center (y-direction)
+
+        :return: x and y coordinates in NN coordinate system
+        """
+        # define rotation matrix to rotate back into alignment
+        counterrotation = -ellipse_pa_to_imagingx_angle
+        cd1_1 = np.cos(counterrotation)
+        cd1_2 = -np.sin(counterrotation)
+        cd2_1 = np.sin(counterrotation)
+        cd2_2 = np.cos(counterrotation)
+        # rotation matrix, applied to matching centers ()
+        rotation_by_ellipse_angle = np.array([[cd1_1, cd1_2], [cd2_1, cd2_2]]) * (deltapix_imaging / deltapix_kin_nn)
+
+        kin_nn_x_at_imagingcenter, kin_nn_y_at_imagingcenter = rotation_by_ellipse_angle.dot(
+            np.array([-npix_imaging / 2 - offsetx / deltapix_imaging, -npix_imaging / 2 - offsety / deltapix_imaging])) + [
+                                                           npix_kin_nn / 2, npix_kin_nn / 2]
+
+        kin_nn_x, kin_nn_y = rotation_by_ellipse_angle.dot(np.array([imaging_x, imaging_y]))
+        return kin_nn_x + kin_nn_x_at_imagingcenter, kin_nn_y + kin_nn_y_at_imagingcenter
+
+    def plot_contour_and_grid(self, xcoords, ycoords, orig_image, color, alpha=0.4):
+        """
+        plotting function for visualization of a grid with a single ellipse contour
+
+        :param xcoords: grid x coordinates, flattened
+        :param ycoords: grid y coordinates, flattened
+        :param orig_image: original 2D image
+        :param color: plotting color
+        :param alpha: plotting alpha
+        """
+        ell_cond_2d = np.isclose(orig_image, 0.1, rtol=5e-02)
+        ell_cond = ell_cond_2d.flatten()
+        plt.scatter(xcoords[ell_cond], ycoords[ell_cond], color=color, alpha=alpha, s=10)
+        plt.scatter(xcoords, ycoords, color=color, alpha=alpha, s=1)
+
+    def spectragrid_in_radec(self):
+        """
+        calculates ra and dec coordinates of spectra input coordinate grid
+
+        :return: ra and dec coordinates
+        """
+        spectra_x, spectra_y = self.pix_coords(self.spectra_data, flatten=True)
+        spectra_ra, spectra_dec = self.xy_to_radec(spectra_x, spectra_y, self.spectra_data['transform_pix2angle'],
+                                             self.spectra_data['ra_at_xy0'], self.spectra_data['dec_at_xy0'])
+        return spectra_ra, spectra_dec
+
+    def spectragrid_in_imagingxy(self):
+        """
+        calculates x and y coordinates in the imaging coordinate system of the original spectra input grid
+
+        :return: x and y coordinates
+        """
+        spectra_ra, spectra_dec = self.spectragrid_in_radec()
+        spectra_coords_in_imaging_x, spectra_coords_in_imaging_y = self.radec_to_xy(spectra_ra, spectra_dec,
+                                                                      self.imaging_data['transform_pix2angle'],
+                                                                      self.imaging_data['ra_at_xy0'],
+                                                                      self.imaging_data['dec_at_xy0'])
+        return spectra_coords_in_imaging_x, spectra_coords_in_imaging_y
+
+    def spectragrid_in_kin_nn_xy(self):
+        """
+        calculates x and y coordinates in the NN coordinate system of the original spectra input grid
+
+        :return: x and y coordinates
+        """
+        spectra_coords_in_imaging_x, spectra_coords_in_imaging_y = self.spectragrid_in_imagingxy()
+        kin_nn_x, kin_nn_y = self.rotate_imaging_into_kin_nn(spectra_coords_in_imaging_x, spectra_coords_in_imaging_y,
+                                                             self.imaging_data['ellipse_PA'], self.imaging_deltapix,
+                                                             self.kinNN_data['deltaPix'], self.imaging_data['npix'],
+                                                             self.kinNN_data['npix'],
+                                                             offsetx=self.imaging_data['offset_x'],
+                                                             offsety=self.imaging_data['offset_y'])
+        return kin_nn_x, kin_nn_y
+
+    def interp_image(self):
+        """
+        interpolates kinNN image at the coordinates of the transformed spectra grid
+
+        :return: interpolated image which lines up with spectra coordinates
+        """
+        spectra_kin_nn_x, spectra_kin_nn_y = self.spectragrid_in_kin_nn_xy()
+        x_axis = np.arange(self.kinNN_data['npix'])
+        y_axis = np.arange(self.kinNN_data['npix'])
+        interp_fcn = RectBivariateSpline(x_axis, y_axis, self.kinNN_data['image'])
+        # y and x are flipped in RectBivariateSpline call:
+        flux_interp = interp_fcn.ev(spectra_kin_nn_y, spectra_kin_nn_x).reshape(self.spectra_data['npix'],
+                                                                          self.spectra_data['npix'])
+        return flux_interp


Test patch:
diff --git a/test/test_Data/test_kinematic_data_2D.py b/test/test_Data/test_kinematic_data_2D.py
new file mode 100644
index 000000000..0eaed5ef6
--- /dev/null
+++ b/test/test_Data/test_kinematic_data_2D.py
@@ -0,0 +1,43 @@
+import pytest
+import numpy as np
+import numpy.testing as npt
+
+import lenstronomy.Util.util as util
+
+from lenstronomy.Data.psf import PSF
+from lenstronomy.Data.kinematic_bin_2D import KinBin
+
+import lenstronomy.Util.kernel_util as kernel_util
+
+
+class TestKinBin(object):
+    def setup(self):
+        self.numPix = 10
+        self.num_bin = 4
+        kwargs_kin = {'bin_data': np.zeros(self.num_bin) + 10.,
+                      'bin_cov': np.diag(np.ones(self.num_bin) * 2. ** 2),
+                      'bin_mask': np.zeros((self.numPix, self.numPix))}
+        kwargs_kin['bin_mask'][0, 0] = 1
+        kwargs_kin['bin_mask'][0, 1] = 1
+        kwargs_kin['bin_mask'][1:5, 2:4] = 2
+        kwargs_kin['bin_mask'][9, 9] = 3
+
+        kwargs_kin['bin_data'][1] = 20.
+        kwargs_kin['transform_pix2angle'] = np.array([[1, 0], [0, 1]])
+        kwargs_kin['ra_at_xy_0'] = 0
+        kwargs_kin['dec_at_xy_0'] = 0
+
+        kernel_point_source = kernel_util.kernel_gaussian(num_pix=9, delta_pix=1., fwhm=2.)
+        kwargs_pixel = {'psf_type': 'PIXEL', 'kernel_point_source': kernel_point_source}
+        self.PSF = PSF(**kwargs_pixel)
+        self.kwargs_kin=kwargs_kin
+        self.KinBin = KinBin(psf_class=self.PSF, **kwargs_kin)
+
+    def test_binned_image(self):
+        assert self.KinBin.binned_image(self.kwargs_kin['bin_data'],self.kwargs_kin['bin_mask'])[0, 0] == 20.
+
+    def test_KinBin2kwargs(self):
+        assert self.KinBin.kin_bin2kwargs()['deltaPix'] == 1.
+
+    def test_kin_grid(self):
+        assert np.shape(self.KinBin.kin_grid()[0]) == (self.numPix, self.numPix)
diff --git a/test/test_LensModel/test_Profiles/test_epl.py b/test/test_LensModel/test_Profiles/test_epl.py
index 2a524c146..43185d2d4 100644
--- a/test/test_LensModel/test_Profiles/test_epl.py
+++ b/test/test_LensModel/test_Profiles/test_epl.py
@@ -193,6 +193,39 @@ def test_epl_pemd_convention(self):
                             npt.assert_almost_equal(f_x, f_x_pemd, decimal=4)
                             npt.assert_almost_equal(f_y, f_y_pemd, decimal=4)
 
+class TestEPLQPhi(object):
+    """
+    tests if EPL and EPL_Q_PHI return same values
+    """
+
+    def setup_method(self):
+        from lenstronomy.LensModel.Profiles.epl import EPL
+        self.epl = EPL()
+        from lenstronomy.LensModel.Profiles.epl import EPLQPhi
+        self.epl_qphi = EPLQPhi()
+
+
+    def test_epl_eplqphi(self):
+
+        x, y = util.make_grid(numPix=10, deltapix=0.2)
+
+        theta_E_list = [0.5, 1, 2]
+        gamma_list = [1.8, 2., 2.2]
+        e1_list = [-0.2, 0., 0.2]
+        e2_list = [-0.2, 0., 0.2]
+        for gamma in gamma_list:
+            for e1 in e1_list:
+                for e2 in e2_list:
+                    phi, q = param_util.ellipticity2phi_q(e1, e2)
+                    for theta_E in theta_E_list:
+                        kwargs = {'theta_E': theta_E, 'gamma': gamma, 'e1': e1, 'e2': e2, 'center_x': 0, 'center_y': 0}
+                        f_x, f_y = self.epl.derivatives(x, y, **kwargs)
+                        kwargs_qphi = {'theta_E': theta_E, 'gamma': gamma, 'q': q, 'phi': phi, 'center_x': 0,
+                                       'center_y': 0}
+                        f_x_qphi, f_y_qphi = self.epl_qphi.derivatives(x, y, **kwargs_qphi)
+
+                        npt.assert_almost_equal(f_x, f_x_qphi, decimal=4)
+                        npt.assert_almost_equal(f_y, f_y_qphi, decimal=4)
 
 if __name__ == '__main__':
     pytest.main()
diff --git a/test/test_LensModel/test_Profiles/test_pemd.py b/test/test_LensModel/test_Profiles/test_pemd.py
index 3bf62fce0..798ee7c5d 100644
--- a/test/test_LensModel/test_Profiles/test_pemd.py
+++ b/test/test_LensModel/test_Profiles/test_pemd.py
@@ -228,3 +228,4 @@ def test_density_lens(self):
 
 if __name__ == '__main__':
     pytest.main()
+
diff --git a/test/test_LensModel/test_numeric_lens_differentials.py b/test/test_LensModel/test_numeric_lens_differentials.py
index d7641a17e..a6f5e177e 100644
--- a/test/test_LensModel/test_numeric_lens_differentials.py
+++ b/test/test_LensModel/test_numeric_lens_differentials.py
@@ -282,6 +282,11 @@ def test_EPL(self):
         lens_model = ['EPL']
         self.assert_differentials(lens_model, kwargs)
 
+    def test_EPL_Q_PHI(self):
+        kwargs = {'theta_E': 2., 'q': 0.9, 'phi': 0.1, 'gamma': 2.13}
+        lens_model = ['EPL_Q_PHI']
+        self.assert_differentials(lens_model, kwargs)
+
     def test_EPL_BOXYDISKY(self):
         kwargs = {'theta_E': 2., 'e1': 0.1, 'e2': 0.2, 'gamma': 2.13, 'a_m': 0.1}
         lens_model = ['EPL_BOXYDISKY']
diff --git a/test/test_LightModel/test_Profiles/test_sersic.py b/test/test_LightModel/test_Profiles/test_sersic.py
index 5102564ba..c97c195f3 100644
--- a/test/test_LightModel/test_Profiles/test_sersic.py
+++ b/test/test_LightModel/test_Profiles/test_sersic.py
@@ -1,7 +1,7 @@
 __author__ = 'sibirrer'
 
 
-from lenstronomy.LightModel.Profiles.sersic import Sersic, SersicElliptic, CoreSersic
+from lenstronomy.LightModel.Profiles.sersic import Sersic, SersicElliptic, SersicElliptic_qPhi, CoreSersic
 import lenstronomy.Util.param_util as param_util
 from lenstronomy.Util import util
 import numpy as np
@@ -16,6 +16,7 @@ class TestSersic(object):
     def setup_method(self):
         self.sersic = Sersic(smoothing=0.02)
         self.sersic_elliptic = SersicElliptic(smoothing=0.02, sersic_major_axis=True)
+        self.sersic_elliptic_qphi = SersicElliptic_qPhi(smoothing=0.02, sersic_major_axis=True)
         self.core_sersic = CoreSersic(smoothing=0.02, sersic_major_axis=True)
 
     def test_sersic(self):
@@ -148,6 +149,20 @@ def test_total_flux(self):
         flux_numeric_ell = np.sum(flux_grid) * deltapix ** 2
         npt.assert_almost_equal(flux_numeric_ell / flux_analytic_ell, 1, decimal=2)
 
+    def test_qphi(self):
+        x = np.array([1])
+        y = np.array([2])
+        I0_sersic = 1
+        R_sersic = 1
+        n_sersic = 1
+        phi_G = 1
+        q = 0.9
+        e1, e2 = param_util.phi_q2_ellipticity(phi_G, q)
+        center_x = 0
+        center_y = 0
+        values = self.sersic_elliptic.function(x, y, I0_sersic, R_sersic, n_sersic, e1, e2, center_x, center_y)
+        values_qphi=self.sersic_elliptic_qphi.function(x, y, I0_sersic, R_sersic, n_sersic, q, phi_G, center_x, center_y)
+        npt.assert_almost_equal(values, values_qphi, decimal=6)
 
 if __name__ == '__main__':
     pytest.main()
diff --git a/test/test_Sampling/test_Likelihoods/test_kinematic_2D_likelihood.py b/test/test_Sampling/test_Likelihoods/test_kinematic_2D_likelihood.py
new file mode 100644
index 000000000..adec8a733
--- /dev/null
+++ b/test/test_Sampling/test_Likelihoods/test_kinematic_2D_likelihood.py
@@ -0,0 +1,191 @@
+import numpy as np
+import numpy.testing as npt
+from lenstronomy.Sampling.Likelihoods.kinematic_2D_likelihood import KinLikelihood
+from lenstronomy.LensModel.lens_model import LensModel
+from lenstronomy.LightModel.light_model import LightModel
+from lenstronomy.Data.psf import PSF
+from lenstronomy.Data.kinematic_bin_2D import KinBin
+from lenstronomy.Data.pixel_grid import PixelGrid
+import lenstronomy.Util.kernel_util as kernel_util
+from lenstronomy.Sampling.Likelihoods import kinematic_NN_call
+
+
+class TestKinLikelihood(object):
+    def setup(self):
+        self.kinematic_NN = kinematic_NN_call.KinematicNN().SKiNN_installed
+        # initialize KinLikelihood class and calculate example vrms for testing
+        major_axis_pa = 30 * np.pi / 180
+        mass_q = 0.7
+        light_q = 0.6
+        center_x = 0.0
+        center_y = 0.0
+        self.kwargs_lens = [
+            {'theta_E': 1.5, 'gamma': 2, 'q': mass_q, 'phi': major_axis_pa, 'center_x': center_x, 'center_y': center_y}]
+        self.kwargs_lens_light = [
+            {'amp': 10, 'R_sersic': 1., 'q': light_q, 'phi': major_axis_pa, 'n_sersic': 3., 'center_x': center_x,
+             'center_y': center_y}]
+        self.kwargs_special = {'D_dt': 4000, 'b_ani': 0.2, 'incli': np.pi / 2, 'D_d': 1000}
+
+        lens_model_list = ['EPL_Q_PHI']
+        self.lensModel = LensModel(lens_model_list=lens_model_list, z_lens=0.5)
+        lens_light_model_list = ['SERSIC_ELLIPSE_Q_PHI']
+
+        self.lensLightModel = LightModel(lens_light_model_list)
+        npix = 3
+        delta_pix_image = 0.2
+        transform_pix2angle = np.array([[1, 0], [0, 1]]) * delta_pix_image
+        ra_at_xy_0 = -(npix - 1) / 2. * delta_pix_image
+        dec_at_xy_0 = -(npix - 1) / 2. * delta_pix_image
+
+        image_pixel_grid = PixelGrid(npix, npix, transform_pix2angle, ra_at_xy_0, dec_at_xy_0)
+        self.image_data = self.lensLightModel.surface_brightness(image_pixel_grid._x_grid, image_pixel_grid._y_grid,
+                                                                 self.kwargs_lens_light)
+
+        binmap = np.array([[0, 0, 1], [0, 0, 1], [2, 2, 3]])
+        binned_dummy_data = np.ones(4)
+        delta_pix_kin = delta_pix_image
+        npix_kin = npix
+
+        # create PSF for kinematics
+        sharp_kin_psf = kernel_util.kernel_gaussian(num_pix=npix, delta_pix=delta_pix_image, fwhm=0.001)
+
+        # KinLikelihood inputs: light image kwargs, kinematic image kwargs, PSF, lens model class, light model class
+        self.kwargs_data = {'image_data': self.image_data,
+                            'noise_map': 100 * np.ones((npix, npix)),
+                            'transform_pix2angle': transform_pix2angle,
+                            'ra_at_xy_0': ra_at_xy_0, 'dec_at_xy_0': dec_at_xy_0}
+        self.kwargs_kin = {'bin_data': binned_dummy_data,
+                           'bin_cov': np.diag((binned_dummy_data * 0.05) ** 2),  # 5% error
+                           'bin_mask': binmap,
+                           'ra_at_xy_0': -(npix_kin - 1) / 2. * delta_pix_kin,
+                           'dec_at_xy_0': -(npix_kin - 1) / 2. * delta_pix_kin,
+                           'transform_pix2angle': np.array([[1, 0], [0, 1]]) * delta_pix_kin}
+
+        kwargs_pixelkin = {'psf_type': 'PIXEL', 'kernel_point_source': sharp_kin_psf}
+        self.kinPSF = PSF(**kwargs_pixelkin)
+        _KinBin = KinBin(psf_class=self.kinPSF, **self.kwargs_kin)
+        self._KinLikelihood = KinLikelihood(_KinBin, self.lensModel, self.lensLightModel, self.kwargs_data, idx_lens=0,
+                                            idx_lens_light=0)
+        if self.kinematic_NN:
+            # if NN enabled, calculate realistic values to use as truth for logL tests
+            self.truth_vrms = self._KinLikelihood.calc_vrms(self.kwargs_lens, self.kwargs_lens_light,
+                                                            self.kwargs_special, verbose=False)
+            self.kwargs_kin = {'bin_data': self.truth_vrms,
+                               'bin_cov': np.diag((self.truth_vrms * 0.05) ** 2),  # 5% error
+                               'bin_mask': binmap,
+                               'ra_at_xy_0': -(npix_kin - 1) / 2. * delta_pix_kin,
+                               'dec_at_xy_0': -(npix_kin - 1) / 2. * delta_pix_kin,
+                               'transform_pix2angle': np.array([[1, 0], [0, 1]]) * delta_pix_kin}
+            _KinBin = KinBin(psf_class=self.kinPSF, **self.kwargs_kin)
+            self._KinLikelihood = KinLikelihood(_KinBin, self.lensModel, self.lensLightModel, self.kwargs_data,
+                                                idx_lens=0, idx_lens_light=0)
+
+    def test_logL(self):
+        if self.kinematic_NN:
+            kwargs_lens_close = [{}]
+            kwargs_lens_far = [{}]
+            for key, val in self.kwargs_lens[0].items():
+                kwargs_lens_close[0][key] = val + 0.01
+                kwargs_lens_far[0][key] = val + 0.2
+            kwargs_lens_light_close = [{}]
+            kwargs_lens_light_far = [{}]
+            for key, val in self.kwargs_lens_light[0].items():
+                kwargs_lens_light_close[0][key] = val + 0.01
+                kwargs_lens_light_far[0][key] = val + 0.2
+
+            close_logL = self._KinLikelihood.logL(kwargs_lens_close, kwargs_lens_light_close, self.kwargs_special,
+                                                  verbose=False)
+            far_logL = self._KinLikelihood.logL(kwargs_lens_far, kwargs_lens_light_far, self.kwargs_special,
+                                                verbose=False)
+            assert close_logL > -10  # instance close to truth should be likely (usually around -0.5)
+            assert close_logL > far_logL  # param instance close to truth is more likely than far away
+
+            # test  that doubling variance decreases logL by factor of 2
+            kwargs_kin = self.kwargs_kin.copy()
+            kwargs_kin['bin_cov'] = self.kwargs_kin['bin_cov'] * 2
+            _KinBin = KinBin(psf_class=self.kinPSF, **kwargs_kin)
+            var_kin_likelihood = KinLikelihood(_KinBin, self.lensModel, self.lensLightModel, self.kwargs_data,
+                                               idx_lens=0,
+                                               idx_lens_light=0)
+            var_logL = var_kin_likelihood.logL(kwargs_lens_close, kwargs_lens_light_close, self.kwargs_special,
+                                               verbose=False)
+            npt.assert_almost_equal(var_logL, close_logL / 2, decimal=2)
+
+            # test that if velocity off by 1 sigma, logL off by 1 sigma
+            kwargs_kin = self.kwargs_kin.copy()
+            kwargs_kin['bin_data'][-1] = kwargs_kin['bin_data'][-1] * 1.05  # increase one bin by 1 sigma
+
+            _KinBin = KinBin(psf_class=self.kinPSF, **kwargs_kin)
+            sig_kin_likelihood = KinLikelihood(_KinBin, self.lensModel, self.lensLightModel, self.kwargs_data,
+                                               idx_lens=0,
+                                               idx_lens_light=0)
+            sig_logL = sig_kin_likelihood.logL(kwargs_lens_close, kwargs_lens_light_close, self.kwargs_special,
+                                               verbose=False)
+            npt.assert_almost_equal(sig_logL, close_logL - 1 / 2, decimal=1)
+            kwargs_out_of_bounds=kwargs_lens_close.copy()
+            kwargs_out_of_bounds[0]['theta_E'] = 10
+            out_of_bounds_logL = self._KinLikelihood.logL(kwargs_lens_close, kwargs_lens_light_close, self.kwargs_special, verbose=False)
+            assert out_of_bounds_logL == -10**8
+
+    def test_convert_to_nn_params(self):
+        kwargs_lens_test = [{'theta_E': 2., 'gamma': 2., 'q': 1., 'phi': 0, 'center_x': 0, 'center_y': 0},
+                            {'gamma1': 0.06, 'gamma2': -0.03}]
+        kwargs_lens_light_test = [{'amp': 10, 'R_sersic': 1., 'q': 1, 'phi': 0, 'n_sersic': 3., 'center_x': 0.,
+                                   'center_y': 0.}]
+        params, same_orientation = self._KinLikelihood.convert_to_nn_params(kwargs_lens_test, kwargs_lens_light_test,
+                                                                            self.kwargs_special)
+        npt.assert_array_equal(params, np.array([1., 1., 2., 3., 1., 8.0e-2, 0.5, self.kwargs_special['b_ani'],
+                                                 self.kwargs_special['incli'] * 180 / np.pi]))
+        # check that wrong orientations are caught
+        assert same_orientation == True
+        kwargs_lens_test = [{'theta_E': 2., 'gamma': 2., 'q': 0.5, 'phi': 1.0, 'center_x': 0, 'center_y': 0},
+                            {'gamma1': 0.06, 'gamma2': -0.03}]
+        params, same_orientation = self._KinLikelihood.convert_to_nn_params(kwargs_lens_test, kwargs_lens_light_test,
+                                                                            self.kwargs_special)
+        assert same_orientation == False
+
+    def test_rescale_distance(self):
+        # scale=fiducial
+        kwargs_special = {'D_dt': 2886.544, 'b_ani': 0.1, 'incli': 0., 'D_d': 1215.739}
+        rescaled_map = self._KinLikelihood.rescale_distance(self.image_data, kwargs_special)
+        npt.assert_allclose(rescaled_map, self.image_data, atol=10 ** -4)
+        # scale D_dt
+        kwargs_special = {'D_dt': 2886.544 * 2, 'b_ani': 0.1, 'incli': 0., 'D_d': 1215.739}
+        rescaled_map = self._KinLikelihood.rescale_distance(self.image_data, kwargs_special)
+        npt.assert_allclose(np.sqrt(2) * self.image_data, rescaled_map, atol=10 ** -4)
+        # scale D_d
+        kwargs_special = {'D_dt': 2886.544, 'b_ani': 0.1, 'incli': 0., 'D_d': 1215.739 * 2}
+        rescaled_map = self._KinLikelihood.rescale_distance(self.image_data, kwargs_special)
+        npt.assert_allclose(1 / np.sqrt(2) * self.image_data, rescaled_map, atol=10 ** -4)
+
+    def test_convert_kwargs_to_kinnalign_input(self):
+        self._KinLikelihood.update_image_input(self.kwargs_lens)  # adds PA and centers to kwargs
+        assert 'ellipse_PA' in self._KinLikelihood.image_input
+
+    def test_auto_binning(self):
+        # light weighting testing, change light image
+        dummy_vrms_map = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
+        sharp_image = np.array([[0, 0, 0], [0, 1, 0], [0, 0, 0]]) + 0.0001
+
+        # auto_binning inherits the psf, bin map, and data shape
+        sharp_kin_psf = kernel_util.kernel_gaussian(num_pix=5, delta_pix=0.2, fwhm=0.01)
+        kwargs_pixelkin = {'psf_type': 'PIXEL', 'kernel_point_source': sharp_kin_psf}
+        kin_psf = PSF(**kwargs_pixelkin)
+        _KinBin = KinBin(psf_class=kin_psf, **self.kwargs_kin)
+        _KinLikelihood = KinLikelihood(_KinBin, self.lensModel, self.lensLightModel, self.kwargs_data, idx_lens=0,
+                                       idx_lens_light=0)
+        vrms = _KinLikelihood.auto_binning(dummy_vrms_map, sharp_image)
+        npt.assert_allclose(vrms, [5, 4.5, 7.5, 9], rtol=1e-3)
+        vrms = _KinLikelihood.auto_binning(dummy_vrms_map, np.ones((3, 3)))
+        npt.assert_allclose(vrms, [3, 4.5, 7.5, 9], rtol=1e-3)
+
+        wide_kin_psf = kernel_util.kernel_gaussian(num_pix=5, delta_pix=0.2, fwhm=100)
+        kwargs_pixelkin = {'psf_type': 'PIXEL', 'kernel_point_source': wide_kin_psf}
+        kin_psf = PSF(**kwargs_pixelkin)
+        _KinBin = KinBin(psf_class=kin_psf, **self.kwargs_kin)
+        _KinLikelihood = KinLikelihood(_KinBin, self.lensModel, self.lensLightModel, self.kwargs_data, idx_lens=0,
+                                       idx_lens_light=0)
+
+        vrms = _KinLikelihood.auto_binning(dummy_vrms_map, sharp_image)
+        expected = np.mean(dummy_vrms_map) / np.mean(sharp_image) * np.ones(4) / 25.
+        npt.assert_allclose(vrms, expected, rtol=1e-3)
diff --git a/test/test_Sampling/test_Likelihoods/test_kinematic_NN_call.py b/test/test_Sampling/test_Likelihoods/test_kinematic_NN_call.py
new file mode 100644
index 000000000..a4e14815e
--- /dev/null
+++ b/test/test_Sampling/test_Likelihoods/test_kinematic_NN_call.py
@@ -0,0 +1,56 @@
+from lenstronomy.Sampling.Likelihoods.kinematic_NN_call import KinematicNN
+import numpy as np
+import numpy.testing as npt
+import copy
+import matplotlib.pyplot as plt
+
+
+class TestKinNNCall(object):
+    def setup(self):
+        self.example_input = np.array([8.56306356e-01, 6.77923305e-01, 1.46690845e+00, 3.00871520e+00,
+                                       1.30421229e+00, 8.0e-2, 4.55349851e-01, -5.68615913e-02,
+                                       8.07435011e+01])
+        self.kinematic_NN = KinematicNN()
+
+    def test_generate_map(self):
+        if self.kinematic_NN.SKiNN_installed:
+            map = self.kinematic_NN.generate_map(input_p=self.example_input)
+            # check that returned map is 551x551 pixels at approximately 200 km/s
+            npt.assert_equal(np.shape(map), (551, 551))
+            npt.assert_allclose(np.mean(map), 235.19, rtol=1e-2)
+            # plt.imshow(map)
+            # plt.colorbar()
+            # plt.show()
+
+    def test_check_bounds(self):
+        if self.kinematic_NN.SKiNN_installed:
+            verbose = True
+            in_bounds = self.kinematic_NN.check_bounds(self.example_input, same_orientation=True, verbose=verbose)
+            assert (in_bounds is True)
+            # check if different orientations
+            in_bounds = self.kinematic_NN.check_bounds(self.example_input, same_orientation=False, verbose=verbose)
+            assert (in_bounds is False)
+            # check all params out of bounds
+            for i in range(len(self.example_input)):
+                test_input = copy.copy(self.example_input)
+                test_input[i] = 100
+                in_bounds = self.kinematic_NN.check_bounds(test_input, same_orientation=True, verbose=verbose)
+                assert (in_bounds is False)
+
+            # check r_eff > theta_E
+            test_input = copy.copy(self.example_input)
+            test_input[4] = 1.6
+            in_bounds = self.kinematic_NN.check_bounds(test_input, same_orientation=True, verbose=verbose)
+            assert (in_bounds is False)
+
+            # check inclination
+            test_input = copy.copy(self.example_input)
+            test_input[8] = 10
+            in_bounds = self.kinematic_NN.check_bounds(test_input, same_orientation=True, verbose=verbose)
+            assert (in_bounds is False)
+
+    def test_plot_map(self):
+        if self.kinematic_NN.SKiNN_installed:
+            f, ax = plt.subplots(1, 1, figsize=(4, 4))
+            ax = self.kinematic_NN.plot_map(ax=ax, input_p=self.example_input)
+            plt.close()
diff --git a/test/test_Sampling/test_likelihood.py b/test/test_Sampling/test_likelihood.py
index 708a1a96a..05f10c730 100644
--- a/test/test_Sampling/test_likelihood.py
+++ b/test/test_Sampling/test_likelihood.py
@@ -10,6 +10,9 @@
 from lenstronomy.Sampling.parameters import Param
 from lenstronomy.Data.imaging_data import ImageData
 from lenstronomy.Data.psf import PSF
+from lenstronomy.Data.kinematic_bin_2D import KinBin
+from lenstronomy.Sampling.Likelihoods import kinematic_NN_call
+import lenstronomy.Util.kernel_util as kernel_util
 
 
 class TestLikelihoodModule(object):
@@ -42,7 +45,7 @@ def setup_method(self):
         kwargs_spep = {'theta_E': 1., 'gamma': 1.95, 'center_x': 0, 'center_y': 0, 'e1': 0.1, 'e2': 0.1}
 
         self.kwargs_lens = [kwargs_spep]
-        kwargs_sersic = {'amp': 1/0.05**2., 'R_sersic': 0.1, 'n_sersic': 2, 'center_x': 0, 'center_y': 0}
+        kwargs_sersic = {'amp': 1 / 0.05 ** 2., 'R_sersic': 0.1, 'n_sersic': 2, 'center_x': 0, 'center_y': 0}
         # 'SERSIC_ELLIPSE': elliptical Sersic profile
         kwargs_sersic_ellipse = {'amp': 1., 'R_sersic': .6, 'n_sersic': 3, 'center_x': 0, 'center_y': 0}
 
@@ -52,11 +55,13 @@ def setup_method(self):
                            'source_amp': 1.}]  # quasar point source position in the source plane and intrinsic brightness
         self.kwargs_cosmo = {'D_dt': 1000}
         kwargs_numerics = {'supersampling_factor': 1, 'supersampling_convolution': False}
-        lens_model_class, source_model_class, lens_light_model_class, point_source_class, extinction_class = class_creator.create_class_instances(**kwargs_model)
+        lens_model_class, source_model_class, lens_light_model_class, point_source_class, extinction_class = class_creator.create_class_instances(
+            **kwargs_model)
         imageModel = ImageModel(data_class, psf_class, lens_model_class, source_model_class,
-                                lens_light_model_class, point_source_class, extinction_class, kwargs_numerics=kwargs_numerics)
+                                lens_light_model_class, point_source_class, extinction_class,
+                                kwargs_numerics=kwargs_numerics)
         image_sim = sim_util.simulate_simple(imageModel, self.kwargs_lens, self.kwargs_source,
-                                         self.kwargs_lens_light, self.kwargs_ps)
+                                             self.kwargs_lens_light, self.kwargs_ps)
         ra_pos, dec_pos = imageModel.PointSource.image_position(kwargs_ps=self.kwargs_ps, kwargs_lens=self.kwargs_lens)
 
         data_class.update_data(image_sim)
@@ -70,15 +75,16 @@ def setup_method(self):
             'supersampling_convolution': False}
 
         kwargs_constraints = {
-                                   'num_point_source_list': [4],
-                                   'solver_type': 'NONE',  # 'PROFILE', 'PROFILE_SHEAR', 'ELLIPSE', 'CENTER'
-                                   'Ddt_sampling': True
-                                   }
+            'num_point_source_list': [4],
+            'solver_type': 'NONE',  # 'PROFILE', 'PROFILE_SHEAR', 'ELLIPSE', 'CENTER'
+            'Ddt_sampling': True
+        }
 
-        def condition_definition(kwargs_lens, kwargs_source, kwargs_lens_light, kwargs_ps=None, kwargs_special=None, kwargs_extinction=None):
+        def condition_definition(kwargs_lens, kwargs_source, kwargs_lens_light, kwargs_ps=None, kwargs_special=None,
+                                 kwargs_extinction=None):
             logL = 0
             if kwargs_lens_light[0]['R_sersic'] > kwargs_source[0]['R_sersic']:
-                logL -= 10**15
+                logL -= 10 ** 15
             return logL
 
         kwargs_likelihood = {'force_no_add_image': True,
@@ -94,7 +100,8 @@ def condition_definition(kwargs_lens, kwargs_source, kwargs_lens_light, kwargs_p
                              'custom_logL_addition': condition_definition,
                              'image_position_likelihood': True
                              }
-        self.kwargs_data = {'multi_band_list': [[kwargs_band, kwargs_psf, kwargs_numerics]], 'multi_band_type': 'single-band',
+        self.kwargs_data = {'multi_band_list': [[kwargs_band, kwargs_psf, kwargs_numerics]],
+                            'multi_band_type': 'single-band',
                             'time_delays_measured': np.ones(3),
                             'time_delays_uncertainties': np.ones(3),
                             'flux_ratios': np.ones(3),
@@ -104,31 +111,112 @@ def condition_definition(kwargs_lens, kwargs_source, kwargs_lens_light, kwargs_p
                             }
         self.param_class = Param(self.kwargs_model, **kwargs_constraints)
         self.imageModel = ImageModel(data_class, psf_class, lens_model_class, source_model_class,
-                                lens_light_model_class,
-                                point_source_class, kwargs_numerics=kwargs_numerics)
-        self.Likelihood = LikelihoodModule(kwargs_data_joint=self.kwargs_data, kwargs_model=kwargs_model, param_class=self.param_class, **kwargs_likelihood)
+                                     lens_light_model_class,
+                                     point_source_class, kwargs_numerics=kwargs_numerics)
+        self.Likelihood = LikelihoodModule(kwargs_data_joint=self.kwargs_data, kwargs_model=kwargs_model,
+                                           param_class=self.param_class, **kwargs_likelihood)
         self.kwargs_band = kwargs_band
         self.kwargs_psf = kwargs_psf
         self.numPix = numPix
 
     def test_logL(self):
         args = self.param_class.kwargs2args(kwargs_lens=self.kwargs_lens, kwargs_source=self.kwargs_source,
-                                            kwargs_lens_light=self.kwargs_lens_light, kwargs_ps=self.kwargs_ps, kwargs_special=self.kwargs_cosmo)
+                                            kwargs_lens_light=self.kwargs_lens_light, kwargs_ps=self.kwargs_ps,
+                                            kwargs_special=self.kwargs_cosmo)
 
         logL = self.Likelihood.logL(args, verbose=True)
         num_data_evaluate = self.Likelihood.num_data
-        npt.assert_almost_equal(logL/num_data_evaluate, -1/2., decimal=1)
+        npt.assert_almost_equal(logL / num_data_evaluate, -1 / 2., decimal=1)
 
     def test_time_delay_likelihood(self):
         kwargs_likelihood = {'time_delay_likelihood': True,
                              }
-        likelihood = LikelihoodModule(kwargs_data_joint=self.kwargs_data, kwargs_model=self.kwargs_model, param_class=self.param_class, **kwargs_likelihood)
+        likelihood = LikelihoodModule(kwargs_data_joint=self.kwargs_data, kwargs_model=self.kwargs_model,
+                                      param_class=self.param_class, **kwargs_likelihood)
         args = self.param_class.kwargs2args(kwargs_lens=self.kwargs_lens, kwargs_source=self.kwargs_source,
-                                            kwargs_lens_light=self.kwargs_lens_light, kwargs_ps=self.kwargs_ps, kwargs_special=self.kwargs_cosmo)
+                                            kwargs_lens_light=self.kwargs_lens_light, kwargs_ps=self.kwargs_ps,
+                                            kwargs_special=self.kwargs_cosmo)
 
         logL = likelihood.logL(args, verbose=True)
         npt.assert_almost_equal(logL, -1328.821179288249, decimal=-1)
 
+    def test_kin_likelihood(self):
+        self.kinematic_NN = kinematic_NN_call.KinematicNN().SKiNN_installed
+        if self.kinematic_NN:
+            # since this doesn't work for the SPEP model, test it with a different lens model:
+            kwargs_eplqphi = {'theta_E': 1., 'gamma': 1.95, 'center_x': 0, 'center_y': 0, 'q': 0.7, 'phi': 0.7}
+            kwargs_lens = [kwargs_eplqphi]
+            kwargs_sersic_ellipse_qphi = {'amp': 1., 'R_sersic': .51, 'n_sersic': 2.5, 'center_x': 0, 'center_y': 0,
+                                          'q': 0.99, 'phi': 0.7}
+            kwargs_lens_light = [kwargs_sersic_ellipse_qphi]
+            kwargs_model = {'lens_model_list': ['EPL_Q_PHI'],
+                            'lens_light_model_list': ['SERSIC_ELLIPSE_Q_PHI'],
+                            'source_light_model_list': ['SERSIC'],
+                            'point_source_model_list': ['SOURCE_POSITION'],
+                            'fixed_magnification_list': [True],
+                            'z_lens': 1}
+            kwargs_constraints = {
+                'num_point_source_list': [4],
+                'solver_type': 'NONE',  # 'PROFILE', 'PROFILE_SHEAR', 'ELLIPSE', 'CENTER'
+                'Ddt_sampling': True
+            }
+            kwargs_likelihood = {'kinematic_2d_likelihood': False, 'kin_lens_idx': 0, 'kin_lens_light_idx': 0,
+                                 'time_delay_likelihood': False}
+            param_class = Param(kwargs_model, **kwargs_constraints)
+            Likelihood = LikelihoodModule(kwargs_data_joint=self.kwargs_data, kwargs_model=kwargs_model,
+                                          param_class=param_class, **kwargs_likelihood)
+            kwargs_cosmo = {'D_dt': 2000, 'D_d': 1000, 'b_ani': 0.2, 'incli': np.pi / 2 - 0.01}
+            args = param_class.kwargs2args(kwargs_lens=kwargs_lens, kwargs_source=self.kwargs_source,
+                                           kwargs_lens_light=kwargs_lens_light, kwargs_ps=self.kwargs_ps,
+                                           kwargs_special=kwargs_cosmo)
+            logL_nokin = Likelihood.logL(args, verbose=True)
+            # logL should be same order as SPEP value, buthas different r_sersic so it will be worse
+            npt.assert_allclose(logL_nokin, -1328.821179288249, rtol=1)
+
+            # Now add kinematic likelihood
+            # for simplicity, set kin image data to same as light data
+            numPix = 50  # cutout pixel size
+            deltaPix = 0.1  # pixel size in arcsec (area per pixel = deltaPix**2)
+
+            binmap = np.zeros_like(self.kwargs_band['image_data'])  # one single bin across whole image
+            binned_dummy_data = np.array([200])
+            delta_pix_kin = deltaPix
+            npix_kin = numPix
+
+            kwargs_kin = {'bin_data': binned_dummy_data,
+                          'bin_cov': np.diag((binned_dummy_data * 0.05) ** 2),  # 5% error
+                          'bin_mask': binmap,
+                          'ra_at_xy_0': -(npix_kin - 1) / 2. * delta_pix_kin,
+                          'dec_at_xy_0': -(npix_kin - 1) / 2. * delta_pix_kin,
+                          'transform_pix2angle': np.array([[1, 0], [0, 1]]) * delta_pix_kin}
+
+            kinkernel_point_source = kernel_util.kernel_gaussian(num_pix=9, delta_pix=0.2, fwhm=1.)
+            kwargs_pixelkin = {'psf_type': 'PIXEL', 'kernel_point_source': kinkernel_point_source}
+            kinPSF = PSF(**kwargs_pixelkin)
+            _KinBin = KinBin(psf_class=kinPSF, **kwargs_kin)
+            kwargs_data_kin = self.kwargs_data.copy()  # add kinematics data to kwargs_data
+            kwargs_data_kin['kinematic_data'] = _KinBin
+
+            # confirm that full likelihood is now lens+kin
+            kwargs_likelihood = {'kinematic_2d_likelihood': True, 'kin_lens_idx': 0, 'kin_lens_light_idx': 0,
+                                 'time_delay_likelihood': False}
+
+            kwargs_constraints['kinematic_sampling'] = True
+            param_class = Param(kwargs_model, **kwargs_constraints)
+            Likelihood = LikelihoodModule(kwargs_data_joint=kwargs_data_kin, kwargs_model=kwargs_model,
+                                          param_class=param_class, **kwargs_likelihood)
+            args = param_class.kwargs2args(kwargs_lens=kwargs_lens, kwargs_source=self.kwargs_source,
+                                           kwargs_lens_light=kwargs_lens_light, kwargs_ps=self.kwargs_ps,
+                                           kwargs_special=kwargs_cosmo)
+
+            logL = Likelihood.logL(args, verbose=True)
+            # With only one bin, the new logL should be worse than without kin by
+            # half the mean chi2 averaged over the whole image
+            image_averaged_chi2 = (
+                        (np.mean(Likelihood.kinematic_2D_likelihood.vrms) - binned_dummy_data) ** 2 / kwargs_kin[
+                    'bin_cov']).flatten()
+            npt.assert_almost_equal(logL_nokin - logL, image_averaged_chi2 / 2)
+
     def test_check_bounds(self):
         penalty, bound_hit = self.Likelihood.check_bounds(args=[0, 1], lowerLimit=[1, 0], upperLimit=[2, 2],
                                                           verbose=True)
@@ -136,11 +224,11 @@ def test_check_bounds(self):
 
     def test_pixelbased_modelling(self):
         ss_source = 2
-        numPix_source = self.numPix*ss_source
+        numPix_source = self.numPix * ss_source
         n_scales = 3
         kwargs_pixelbased = {
             'source_interpolation': 'nearest',
-            'supersampling_factor_source': ss_source, # supersampling of pixelated source grid
+            'supersampling_factor_source': ss_source,  # supersampling of pixelated source grid
 
             # following choices are to minimize pixel solver runtime (not to get accurate reconstruction!)
             'threshold_decrease_type': 'none',
@@ -150,9 +238,9 @@ def test_pixelbased_modelling(self):
             'num_iter_weights': 2,
         }
         kwargs_likelihood = {
-             'image_likelihood': True,
-             'kwargs_pixelbased': kwargs_pixelbased,
-             'check_positive_flux': True,  # effectively not applied, activated for code coverage purposes
+            'image_likelihood': True,
+            'kwargs_pixelbased': kwargs_pixelbased,
+            'check_positive_flux': True,  # effectively not applied, activated for code coverage purposes
         }
         kernel = PSF(**self.kwargs_psf).kernel_point_source
         kwargs_psf = {'psf_type': 'PIXEL', 'kernel_point_source': kernel}
@@ -163,27 +251,30 @@ def test_pixelbased_modelling(self):
             'lens_light_model_list': ['SLIT_STARLETS'],
             'source_light_model_list': ['SLIT_STARLETS'],
         }
-        
-        kwargs_fixed_source = [{'n_scales': n_scales, 'n_pixels': numPix_source**2, 'scale': 1, 'center_x': 0, 'center_y': 0}]
-        kwargs_fixed_lens_light = [{'n_scales': n_scales, 'n_pixels': self.numPix**2, 'scale': 1, 'center_x': 0, 'center_y': 0}]
+
+        kwargs_fixed_source = [
+            {'n_scales': n_scales, 'n_pixels': numPix_source ** 2, 'scale': 1, 'center_x': 0, 'center_y': 0}]
+        kwargs_fixed_lens_light = [
+            {'n_scales': n_scales, 'n_pixels': self.numPix ** 2, 'scale': 1, 'center_x': 0, 'center_y': 0}]
         kwargs_constraints = {'source_grid_offset': True}
-        param_class = Param(kwargs_model, 
-                            kwargs_fixed_source=kwargs_fixed_source, 
+        param_class = Param(kwargs_model,
+                            kwargs_fixed_source=kwargs_fixed_source,
                             kwargs_fixed_lens_light=kwargs_fixed_lens_light,
                             **kwargs_constraints)
 
-        likelihood = LikelihoodModule(kwargs_data_joint=kwargs_data, kwargs_model=kwargs_model, 
+        likelihood = LikelihoodModule(kwargs_data_joint=kwargs_data, kwargs_model=kwargs_model,
                                       param_class=param_class, **kwargs_likelihood)
 
-        kwargs_source = [{'amp': np.ones(n_scales*numPix_source**2)}]
-        kwargs_lens_light = [{'amp': np.ones(n_scales*self.numPix**2)}]
+        kwargs_source = [{'amp': np.ones(n_scales * numPix_source ** 2)}]
+        kwargs_lens_light = [{'amp': np.ones(n_scales * self.numPix ** 2)}]
         kwargs_special = {'delta_x_source_grid': 0, 'delta_y_source_grid': 0}
         args = param_class.kwargs2args(kwargs_lens=self.kwargs_lens, kwargs_source=kwargs_source,
                                        kwargs_lens_light=kwargs_lens_light, kwargs_special=kwargs_special)
 
         logL = likelihood.logL(args, verbose=True)
         num_data_evaluate = likelihood.num_data
-        npt.assert_almost_equal(logL/num_data_evaluate, -1/2., decimal=1)
+        npt.assert_almost_equal(logL / num_data_evaluate, -1 / 2., decimal=1)
+
 
 if __name__ == '__main__':
     pytest.main()
diff --git a/test/test_Sampling/test_special_param.py b/test/test_Sampling/test_special_param.py
index 5919ab598..486ab1599 100644
--- a/test/test_Sampling/test_special_param.py
+++ b/test/test_Sampling/test_special_param.py
@@ -11,12 +11,12 @@ class TestParam(object):
     def setup_method(self):
         self.param = SpecialParam(Ddt_sampling=True, kwargs_fixed=None, point_source_offset=True, num_images=2,
                                   source_size=True, num_tau0=2, num_z_sampling=3, source_grid_offset=True,
-                                  kwargs_lower={'z_sampling': [0.05, 0.1, 0.5]},
-                                  kwargs_upper={'z_sampling': [0.2, 1., 1.]}
+				  kinematic_sampling=True, 
+				  kwargs_lower={'z_sampling': [0.05, 0.1, 0.5]}, kwargs_upper={'z_sampling': [0.2, 1., 1.]}
                                   )
         self.kwargs = {'D_dt': 1988, 'delta_x_image': [0, 0], 'delta_y_image': [0, 0], 'source_size': 0.1,
                        'tau0_list': [0, 1], 'z_sampling': np.array([0.1, 0.5, 2]),
-                       'delta_x_source_grid': 0, 'delta_y_source_grid': 0}
+                       'delta_x_source_grid': 0, 'delta_y_source_grid': 0, 'b_ani':0.1, 'incli':0.,'D_d':2000}
 
     def test_get_setParams(self):
         args = self.param.set_params(self.kwargs)
@@ -26,9 +26,12 @@ def test_get_setParams(self):
             npt.assert_almost_equal(args[k], args_new[k], decimal=8)
 
         param_fixed = SpecialParam(Ddt_sampling=True, kwargs_fixed=self.kwargs, point_source_offset=True, num_images=2,
-                                   source_size=True, num_z_sampling=3, num_tau0=2)
+                                   source_size=True, num_z_sampling=3, num_tau0=2,kinematic_sampling=True)
         kwargs_new, i = param_fixed.get_params(args=[], i=0)
         kwargs_new['D_dt'] = self.kwargs['D_dt']
+        kwargs_new['b_ani'] = self.kwargs['b_ani']
+        kwargs_new['incli'] = self.kwargs['incli']
+        kwargs_new['D_d'] = self.kwargs['D_d']
 
         special_param = SpecialParam(num_z_sampling=1, kwargs_lower={'z_sampling': [0.1]},
                                      kwargs_upper={'z_sampling': [0.2]})
@@ -38,7 +41,7 @@ def test_get_setParams(self):
 
     def test_num_params(self):
         num, list = self.param.num_param()
-        assert num == 13
+        assert num == 16
 
     def test_mass_scaling(self):
         kwargs_fixed = {}
diff --git a/test/test_Util/test_image_util.py b/test/test_Util/test_image_util.py
index a5fa6d221..562dd0967 100644
--- a/test/test_Util/test_image_util.py
+++ b/test/test_Util/test_image_util.py
@@ -179,6 +179,29 @@ def test_rotateImage():
     npt.assert_almost_equal(im_rot[2, 2], 1., decimal=10)
     npt.assert_almost_equal(im_rot[2, 1], 0.23931518624017073, decimal=10)
 
+def test_shift_image():
+    img = np.zeros((5, 5))
+    img[2, 2] = 1
+    img[1, 2] = 0.5
+
+    shift = [1, 1]
+    im_shift = image_util.shift_image(img, shift)
+    npt.assert_almost_equal(im_shift[2, 3], 0.5, decimal=10)
+    npt.assert_almost_equal(im_shift[3, 3], 1., decimal=10)
+    npt.assert_almost_equal(im_shift[2, 1], 0., decimal=10)
+
+    shift = [1, 0]
+    im_shift = image_util.shift_image(img, shift)
+    npt.assert_almost_equal(im_shift[2, 2], 0.5, decimal=10)
+    npt.assert_almost_equal(im_shift[3, 2], 1., decimal=10)
+    npt.assert_almost_equal(im_shift[2, 1], 0., decimal=10)
+
+    img = np.zeros((5, 5))
+    img[2,2] = 1
+    shift = [0.5,0.5]
+    im_shift = image_util.shift_image(img, shift)
+    npt.assert_almost_equal(np.sum(im_shift), np.sum(img), decimal=10)
+
 
 def test_re_size_array():
     numPix = 9
diff --git a/test/test_Util/test_kin_sampling_util.py b/test/test_Util/test_kin_sampling_util.py
new file mode 100644
index 000000000..0bab48724
--- /dev/null
+++ b/test/test_Util/test_kin_sampling_util.py
@@ -0,0 +1,83 @@
+import lenstronomy.Util.kin_sampling_util as kin_sampling_util
+
+import numpy as np
+import pytest
+import numpy.testing as npt
+import matplotlib.pyplot as plt
+
+
+class TestKinSamplingUtil(object):
+    def setup(self):
+        self.image = np.array([[0, 1, 0, 0, 0], [1, 1, 1, 1, 1], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0]])
+
+        transform_pix2angle = np.array([[1, 0], [0, 1]])  # linear translation matrix of no rotation
+        ra_at_xy_0, dec_at_xy_0 = transform_pix2angle.dot([-5 / 2, -5 / 2]).T
+        self.kwargs_pixel_rot_imaging = {'nx': 5, 'ny': 5,  # number of pixels per axis
+                                     'ra_at_xy_0': ra_at_xy_0,  # RA at pixel (0,0)
+                                     'dec_at_xy_0': dec_at_xy_0,  # DEC at pixel (0,0)
+                                     'transform_pix2angle': transform_pix2angle}
+        self.imaging_inputs = {'image': self.image, 'transform_pix2angle': self.kwargs_pixel_rot_imaging['transform_pix2angle'],
+                           'ra_at_xy0': self.kwargs_pixel_rot_imaging['ra_at_xy_0'],
+                           'dec_at_xy0': self.kwargs_pixel_rot_imaging['dec_at_xy_0'],
+                           'ellipse_PA': 0,
+                           'offset_x': 0, 'offset_y': 0}
+        self.kinNN_inputs = {'image': self.image, 'deltaPix': 1}
+
+    def define_spectra_kwargs_from_rotation_angle(self, rotation_angle):
+        transform_pix2angle = np.array([[np.cos(rotation_angle), -np.sin(rotation_angle)], [np.sin(rotation_angle),
+                                                                                            np.cos(
+                                                                                                rotation_angle)]])  # linear translation matrix of a 90 degree rotation
+        ra_at_xy_0, dec_at_xy_0 = transform_pix2angle.dot([-5 / 2, -5 / 2]).T
+        kwargs_pixel_rot_spectra = {'nx': 5, 'ny': 5,  # number of pixels per axis
+                                 'ra_at_xy_0': ra_at_xy_0,  # RA at pixel (0,0)
+                                 'dec_at_xy_0': dec_at_xy_0,  # DEC at pixel (0,0)
+                                 'transform_pix2angle': transform_pix2angle}
+        spectra_inputs = {'image': self.image, 'transform_pix2angle': kwargs_pixel_rot_spectra['transform_pix2angle'],
+                       'ra_at_xy0': kwargs_pixel_rot_spectra['ra_at_xy_0'],
+                       'dec_at_xy0': kwargs_pixel_rot_spectra['dec_at_xy_0']}
+        return spectra_inputs
+
+    def test_rotation(self):
+        # spectra grid not rotated from imaging:
+        spectra_inputs = self.define_spectra_kwargs_from_rotation_angle(0)
+        skinn_align = kin_sampling_util.KinNNImageAlign(spectra_inputs, self.imaging_inputs, self.kinNN_inputs)
+        spectra_img_coords = skinn_align.spectragrid_in_imagingxy()
+        imaging_img_coords = skinn_align.pix_coords(self.imaging_inputs)
+        npt.assert_allclose(imaging_img_coords, spectra_img_coords, atol=1e-7)
+
+        spectra_inputs = self.define_spectra_kwargs_from_rotation_angle(2 * np.pi)
+        skinn_align = kin_sampling_util.KinNNImageAlign(spectra_inputs, self.imaging_inputs, self.kinNN_inputs)
+        spectra_img_coords = skinn_align.spectragrid_in_imagingxy()
+        imaging_img_coords = skinn_align.pix_coords(self.imaging_inputs)
+        npt.assert_allclose(imaging_img_coords, spectra_img_coords, atol=1e-7)
+
+        # spectra_grid 90 degrees from imaging:
+        spectra_inputs = self.define_spectra_kwargs_from_rotation_angle(np.pi / 2)
+        skinn_align = kin_sampling_util.KinNNImageAlign(spectra_inputs, self.imaging_inputs, self.kinNN_inputs)
+        spectra_img_coords = skinn_align.spectragrid_in_imagingxy()
+        imaging_img_coords = skinn_align.pix_coords(self.imaging_inputs)
+        npt.assert_allclose(imaging_img_coords[0], spectra_img_coords[1], atol=1e-7)  # old x coordinate is new y coordinate
+
+    def test_samegrid(self):
+        # check interpolation in case where grids are the same
+        spectra_inputs = self.define_spectra_kwargs_from_rotation_angle(0)
+        skinn_align = kin_sampling_util.KinNNImageAlign(spectra_inputs, self.imaging_inputs, self.kinNN_inputs)
+        interp_image = skinn_align.interp_image()
+        npt.assert_almost_equal(interp_image, self.image)
+
+    def test_offset(self):
+        imaging_inputs = {'image': self.image, 'transform_pix2angle': self.kwargs_pixel_rot_imaging['transform_pix2angle'],
+                      'ra_at_xy0': self.kwargs_pixel_rot_imaging['ra_at_xy_0'],
+                      'dec_at_xy0': self.kwargs_pixel_rot_imaging['dec_at_xy_0'],
+                      'ellipse_PA': 0,
+                      'offset_x': 1, 'offset_y': 1}
+        spectra_inputs = self.define_spectra_kwargs_from_rotation_angle(0)
+        skinn_align = kin_sampling_util.KinNNImageAlign(spectra_inputs, imaging_inputs, self.kinNN_inputs)
+        interp_image = skinn_align.interp_image()
+        # #uncomment to see plots
+        # plt.imshow(image)
+        # plt.show()
+        # plt.imshow(interp_image)
+        # plt.show()
+        expected_img = np.array([[0, 0, 1, 0, 0], [0, 0, 1, 0, 0], [1, 1, 1, 1, 1], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0]])
+        npt.assert_allclose(interp_image, expected_img, atol=1e-7)
diff --git a/test_requirements.txt b/test_requirements.txt
index 9f93a5910..aa564fe07 100644
--- a/test_requirements.txt
+++ b/test_requirements.txt
@@ -1,3 +1,4 @@
 # requirements for tests
 colossus
-git+https://github.com/aymgal/coolest@v0.1.0
+git+https://github.com/mattgomer/SKiNN@main#egg=SKiNN
+git+https://github.com/aymgal/coolest@v0.1.0
\ No newline at end of file