How to use the lenstronomy.Util.param_util function in lenstronomy

f_x_e, f_y_e = self.nfw_e.derivatives(x, y, Rs, alpha_Rs, e1, e2)
        npt.assert_almost_equal(f_x[0], f_x_e[0], decimal=5)
        npt.assert_almost_equal(f_y[0], f_y_e[0], decimal=5)
        x = np.array([0])
        y = np.array([0])
        alpha_Rs = 0
        f_x, f_y = self.nfw_e.derivatives(x, y, Rs, alpha_Rs, e1, e2)
        npt.assert_almost_equal(f_x[0], 0, decimal=5)
        npt.assert_almost_equal(f_y[0], 0, decimal=5)

        x = np.array([1,3,4])
        y = np.array([2,1,1])
        alpha_Rs = 1.
        q = .8
        phi_G = 0
        e1, e2 = param_util.phi_q2_ellipticity(phi_G, q)
        values = self.nfw_e.derivatives(x, y, Rs, alpha_Rs, e1, e2)
        npt.assert_almost_equal(values[0][0], 0.32458737284934414, decimal=5)
        npt.assert_almost_equal(values[1][0], 0.9737621185480323, decimal=5)
        npt.assert_almost_equal(values[0][1], 0.76249351329615234, decimal=5)
        npt.assert_almost_equal(values[1][1], 0.38124675664807617, decimal=5)

# 'EXERNAL_SHEAR': external shear
        kwargs_shear = {'gamma1': 0.01, 'gamma2': 0.01}  # gamma_ext: shear strength, psi_ext: shear angel (in radian)
        phi, q = 0.2, 0.8
        e1, e2 = param_util.phi_q2_ellipticity(phi, q)
        kwargs_spemd = {'theta_E': 1., 'gamma': 1.8, 'center_x': 0, 'center_y': 0, 'e1': e1, 'e2': e2}

        lens_model_list = ['SPEP', 'SHEAR']
        self.kwargs_lens = [kwargs_spemd, kwargs_shear]
        lens_model_class = LensModel(lens_model_list=lens_model_list)
        # list of light profiles (for lens and source)
        # 'SERSIC': spherical Sersic profile
        kwargs_sersic = {'amp': 1., 'R_sersic': 0.1, 'n_sersic': 2, 'center_x': 0, 'center_y': 0}
        # 'SERSIC_ELLIPSE': elliptical Sersic profile
        phi, q = 0.2, 0.9
        e1, e2 = param_util.phi_q2_ellipticity(phi, q)
        kwargs_sersic_ellipse = {'amp': 1., 'R_sersic': .6, 'n_sersic': 7, 'center_x': 0, 'center_y': 0,
                                 'e1': e1, 'e2': e2}

        lens_light_model_list = ['SERSIC']
        self.kwargs_lens_light = [kwargs_sersic]
        lens_light_model_class = LightModel(light_model_list=lens_light_model_list)
        source_model_list = ['SERSIC_ELLIPSE']
        self.kwargs_source = [kwargs_sersic_ellipse]
        source_model_class = LightModel(light_model_list=source_model_list)
        self.kwargs_ps = [{'ra_source': 0.0, 'dec_source': 0.0,
                           'source_amp': 10.}]  # quasar point source position in the source plane and intrinsic brightness
        point_source_class = PointSource(point_source_type_list=['SOURCE_POSITION'], fixed_magnification_list=[True])

        kwargs_numerics = {'supersampling_factor': 3, 'supersampling_convolution': False, 'compute_mode': 'regular',
                           'point_source_supersampling_factor': 3}
        imageModel = ImageModel(data_class, psf_class, lens_model_class, source_model_class,

def test_critical_curves(self):
        lens_model_list = ['SPEP']
        phi, q = 1., 0.8
        e1, e2 = param_util.phi_q2_ellipticity(phi, q)
        kwargs_lens = [{'theta_E': 1., 'gamma': 2., 'e1': e1, 'e2': e2, 'center_x': 0, 'center_y': 0}]
        lensModel = LensModelExtensions(LensModel(lens_model_list))
        ra_crit_list, dec_crit_list, ra_caustic_list, dec_caustic_list = lensModel.critical_curve_caustics(kwargs_lens,
                                                                                compute_window=5, grid_scale=0.005)
        print(ra_caustic_list)
        npt.assert_almost_equal(ra_caustic_list[0][3], -0.25629009803139047, decimal=5)
        npt.assert_almost_equal(dec_caustic_list[0][3], -0.39153358367275115, decimal=5)
        npt.assert_almost_equal(ra_crit_list[0][3], -0.53249999999999997, decimal=5)
        npt.assert_almost_equal(dec_crit_list[0][3], -1.2536936868024853, decimal=5)

x = np.array([1])
        y = np.array([2])
        theta_E = 1.
        q = 0.99999
        phi_G = 0
        s = 0.0000001
        e1, e2 = param_util.phi_q2_ellipticity(phi_G, q)
        f_x, f_y = self.nie.derivatives(x, y, theta_E, e1, e2, s_scale=s)
        f_x_spemd, f_y_spemd = self.sis.derivatives(x, y, theta_E)
        npt.assert_almost_equal(f_x, f_x_spemd, decimal=4)
        npt.assert_almost_equal(f_y, f_y_spemd, decimal=4)
        if bool_test is True:
            q = 0.99
            s = 0.000001
            phi_G = 0
            e1, e2 = param_util.phi_q2_ellipticity(phi_G, q)
            f_x, f_y = self.nie.derivatives(x, y, theta_E, e1, e2, s_scale=s)
            gamma = 2.
            f_x_spemd, f_y_spemd = self.spemd.derivatives(x, y, theta_E, gamma, e1, e2, s_scale=s)
            print(f_x/f_x_spemd, 'ratio deflections')
            print(1+(1-q)/2)
            npt.assert_almost_equal(f_x, f_x_spemd, decimal=2)
            npt.assert_almost_equal(f_y, f_y_spemd, decimal=2)

# 'EXERNAL_SHEAR': external shear
        kwargs_shear = {'gamma1': 0.01, 'gamma2': 0.01}  # gamma_ext: shear strength, psi_ext: shear angel (in radian)
        phi, q = 0.2, 0.8
        e1, e2 = param_util.phi_q2_ellipticity(phi, q)
        kwargs_spemd = {'theta_E': 1., 'gamma': 1.8, 'center_x': 0, 'center_y': 0, 'e1': e1, 'e2': e2}

        lens_model_list = ['SPEP', 'SHEAR']
        self.kwargs_lens = [kwargs_spemd, kwargs_shear]
        lens_model_class = LensModel(lens_model_list=lens_model_list)
        # list of light profiles (for lens and source)
        # 'SERSIC': spherical Sersic profile
        kwargs_sersic = {'amp': 1., 'R_sersic': 0.1, 'n_sersic': 2, 'center_x': 0, 'center_y': 0}
        # 'SERSIC_ELLIPSE': elliptical Sersic profile
        phi, q = 0.2, 0.9
        e1, e2 = param_util.phi_q2_ellipticity(phi, q)
        kwargs_sersic_ellipse = {'amp': 1., 'R_sersic': .6, 'n_sersic': 7, 'center_x': 0, 'center_y': 0,
                                 'e1': e1, 'e2': e2}

        lens_light_model_list = ['SERSIC']
        self.kwargs_lens_light = [kwargs_sersic]
        lens_light_model_class = LightModel(light_model_list=lens_light_model_list)
        source_model_list = ['SERSIC_ELLIPSE']
        self.kwargs_source = [kwargs_sersic_ellipse]
        source_model_class = LightModel(light_model_list=source_model_list)
        self.kwargs_ps = [{'ra_source': 0.01, 'dec_source': 0.0,
                       'source_amp': 1.}]  # quasar point source position in the source plane and intrinsic brightness
        point_source_class = PointSource(point_source_type_list=['SOURCE_POSITION'], fixed_magnification_list=[True])
        kwargs_numerics = {'supersampling_factor': 2, 'supersampling_convolution': False}
        imageModel = ImageModel(data_class, psf_class, lens_model_class, source_model_class, lens_light_model_class, point_source_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)

lensModel = LensModel(['SPEP'])
        solver_nfw_ellipse = Solver2Point(lensModel, solver_type='ELLIPSE')
        solver_nfw_center = Solver2Point(lensModel, solver_type='CENTER')
        spep = LensModel(['SPEP'])

        image_position_nfw = LensEquationSolver(LensModel(['SPEP', 'NFW']))
        sourcePos_x = 0.1
        sourcePos_y = 0.03
        deltapix = 0.05
        numPix = 100
        gamma = 1.9
        Rs = 0.1
        nfw = NFW()
        alpha_Rs = nfw._rho02alpha(1., Rs)
        phi_G, q = 0.5, 0.8
        e1, e2 = param_util.phi_q2_ellipticity(phi_G, q)
        kwargs_lens = [{'theta_E': 1., 'gamma': gamma, 'e1': e1, 'e2': e2, 'center_x': 0.1, 'center_y': -0.1},
                       {'Rs': Rs, 'alpha_Rs': alpha_Rs, 'center_x': -0.5, 'center_y': 0.5}]
        x_pos, y_pos = image_position_nfw.findBrightImage(sourcePos_x, sourcePos_y, kwargs_lens, numImages=2, min_distance=deltapix, search_window=numPix*deltapix)
        print(len(x_pos), 'number of images')
        x_pos = x_pos[:2]
        y_pos = y_pos[:2]

        kwargs_init = [{'theta_E': 1, 'gamma': gamma, 'e1': e1, 'e2': e2, 'center_x': 0., 'center_y': 0},
                       {'Rs': Rs, 'alpha_Rs': alpha_Rs, 'center_x': -0.5, 'center_y': 0.5}]
        kwargs_out_center, precision = solver_nfw_center.constraint_lensmodel(x_pos, y_pos, kwargs_init)
        source_x, source_y = spep.ray_shooting(x_pos[0], y_pos[0], kwargs_out_center)
        x_pos_new, y_pos_new = image_position_nfw.findBrightImage(source_x, source_y, kwargs_out_center, numImages=2, min_distance=deltapix, search_window=numPix*deltapix)
        print(kwargs_out_center, 'kwargs_out_center')
        npt.assert_almost_equal(x_pos_new[0], x_pos[0], decimal=2)
        npt.assert_almost_equal(y_pos_new[0], y_pos[0], decimal=2)

def test_function(self):
        """

        :return:
        """
        chameleon = Chameleon()
        triplechameleon = TripleChameleon()

        x = np.linspace(0.1, 10, 10)
        phi_G, q = 0.3, 0.8
        ratio12 = 2.
        ratio13 = 3
        e1, e2 = param_util.phi_q2_ellipticity(phi_G, q)
        kwargs_light = {'amp': 1., 'ratio12': ratio12, 'ratio13': ratio13, 'w_c1': .5, 'w_t1': 1., 'e11': e1, 'e21': e2,
                        'w_c2': .1, 'w_t2': .5, 'e12': e1, 'e22': e2,
                        'w_c3': .1, 'w_t3': .5, 'e13': e1, 'e23': e2
                        }

        amp1 = 1. / (1. + 1. / ratio12 + 1. / ratio13)
        amp2 = amp1 / ratio12
        amp3 = amp1 / ratio13
        kwargs_1 = {'amp': amp1, 'w_c': .5, 'w_t': 1., 'e1': e1, 'e2': e2}
        kwargs_2 = {'amp': amp2, 'w_c': .1, 'w_t': .5, 'e1': e1, 'e2': e2}
        kwargs_3 = {'amp': amp3, 'w_c': .1, 'w_t': .5, 'e1': e1, 'e2': e2}
        flux = triplechameleon.function(x=x, y=1., **kwargs_light)
        flux1 = chameleon.function(x=x, y=1., **kwargs_1)
        flux2 = chameleon.function(x=x, y=1., **kwargs_2)
        flux3 = chameleon.function(x=x, y=1., **kwargs_3)
        npt.assert_almost_equal(flux, flux1 + flux2 + flux3, decimal=8)

:param amp: Amplitude of Gaussian, convention: :math:`A/(2 \pi\sigma^2) \exp(-(x^2+y^2/q^2)/2\sigma^2)`
        :type amp: ``float``
        :param sigma: Standard deviation of Gaussian
        :type sigma: ``float``
        :param e1: Ellipticity parameter 1
        :type e1: ``float``
        :param e2: Ellipticity parameter 2
        :type e2: ``float``
        :param center_x: x coordinate of centroid
        :type center_x: ``float``
        :param center_y: y coordianate of centroid
        :type center_y: ``float``
        :return: Hessian :math:`A/(2 \pi \sigma^2) \exp(-(x^2+y^2/q^2)/2\sigma^2)` for elliptical Gaussian convergence.
        :rtype: tuple ``(float, float, float)`` , or ``(numpy.array, numpy.array, numpy.array)`` with each ``numpy.array``'s shape equal to ``x.shape``.
        """
        phi_g, q = param_util.ellipticity2phi_q(e1, e2)

        if q > 1 - self.min_ellipticity:
            return self.spherical.hessian(x, y, amp, sigma, center_x, center_y)

        # adjusting amplitude to make the notation compatible with the
        # formulae given in Shajib (2019).
        amp_ = amp / (2 * np.pi * sigma**2)

        # converting ellipticity definition from x^2 + y^2/q^2 to q^2*x^2 + y^2
        sigma_ = sigma * q

        x_shift = x - center_x
        y_shift = y - center_y
        cos_phi = np.cos(phi_g)
        sin_phi = np.sin(phi_g)

"""

        :param x: ra-coordinate
        :param y: dec-coordinate
        :param alpha_1: deflection angle at 1 (arcseconds) from the center
        :param w_c: see Suyu+2014
        :param w_t: see Suyu+2014
        :param e1: ellipticity parameter
        :param e2: ellipticity parameter
        :param center_x: ra center
        :param center_y: dec center
        :return: lensing potential
        """

        theta_E_conv, w_c, w_t = self._theta_convert(alpha_1, w_c, w_t)
        phi_G, q = param_util.ellipticity2phi_q(e1, e2)
        s_scale_1 = np.sqrt(4 * w_c ** 2 / (1. + q) ** 2)
        s_scale_2 = np.sqrt(4 * w_t ** 2 / (1. + q) ** 2)
        f_1 = self.nie.function(x, y, theta_E_conv, e1, e2, s_scale_1, center_x, center_y)
        f_2 = self.nie.function(x, y, theta_E_conv, e1, e2, s_scale_2, center_x, center_y)
        f_ = f_1 - f_2
        return f_

x = [e1, e2]
        elif self._solver_type == 'SHAPELETS':
            coeffs = list(kwargs_list[0]['coeffs'])
            [c10, c01] = coeffs[1: 3]
            x = [c10, c01]
        elif self._solver_type == 'THETA_E_PHI':
            theta_E = kwargs_list[0]['theta_E']
            gamma1 = kwargs_list[1]['gamma1']
            gamma2 = kwargs_list[1]['gamma2']
            phi_ext, gamma_ext = param_util.shear_cartesian2polar(gamma1, gamma2)
            x = [theta_E, phi_ext]
        elif self._solver_type == 'THETA_E_ELLIPSE':
            theta_E = kwargs_list[0]['theta_E']
            e1 = kwargs_list[0]['e1']
            e2 = kwargs_list[0]['e2']
            phi_ext, gamma_ext = param_util.shear_cartesian2polar(e1, e2)
            x = [theta_E, phi_ext]
        else:
            raise ValueError("Solver type %s not supported for 2-point solver!" % self._solver_type)
        return x