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

def test_convergence2surface_brightness(self):
        from lenstronomy.LightModel.Profiles.nie import NIE as NIE_Light
        nie_light = NIE_Light()
        kwargs = {'e1': 0.3, 'e2': -0.05, 's_scale': 0.5}
        x, y = util.make_grid(numPix=10, deltapix=0.1)
        f_xx, f_yy, f_xy = self.nie.hessian(x, y, theta_E=1, **kwargs)
        kappa = 1/2. * (f_xx + f_yy)
        flux = nie_light.function(x, y, amp=1, **kwargs)
        npt.assert_almost_equal(kappa/np.sum(kappa), flux/np.sum(flux), decimal=5)

numPix = 100
        deltaPix = 0.05
        kwargs_options = {
                          'lens_light_model_list': ['SERSIC_ELLIPSE', 'SERSIC']}
        kwargs_lens_light = [{'R_sersic': 0.5, 'n_sersic': 4, 'amp': 2, 'e1': 0, 'e2': 0.05},
                             {'R_sersic': 1.5, 'n_sersic': 1, 'amp': 2}]
        lensAnalysis = LensAnalysis(kwargs_options)
        kwargs_interpol = lensAnalysis.light2mass_interpol(lens_light_model_list=['SERSIC_ELLIPSE', 'SERSIC'],
                                                                                          kwargs_lens_light=kwargs_lens_light, numPix=numPix, deltaPix=deltaPix, subgrid_res=1)
        from lenstronomy.LensModel.lens_model import LensModel
        lensModel = LensModel(lens_model_list=['INTERPOL_SCALED'])
        kwargs_lens = [kwargs_interpol]
        import lenstronomy.Util.util as util
        x_grid, y_grid = util.make_grid(numPix, deltapix=deltaPix)
        kappa = lensModel.kappa(x_grid, y_grid, kwargs=kwargs_lens)
        kappa = util.array2image(kappa)
        kappa /= np.mean(kappa)
        flux = lensAnalysis.LensLightModel.surface_brightness(x_grid, y_grid, kwargs_lens_light)
        flux = util.array2image(flux)
        flux /= np.mean(flux)
        #import matplotlib.pyplot as plt
        #plt.matshow(flux-kappa)
        #plt.colorbar()
        #plt.show()
        delta_kappa = (kappa - flux)/flux
        max_delta = np.max(np.abs(delta_kappa))
        assert max_delta < 1
        #assert max_diff < 0.01
        npt.assert_almost_equal(flux[0, 0], kappa[0, 0], decimal=2)

def setup(self):

        # data specifics
        sigma_bkg = 0.01  # background noise per pixel
        exp_time = 100  # exposure time (arbitrary units, flux per pixel is in units #photons/exp_time unit)
        numPix = 100  # cutout pixel size
        deltaPix = 0.05  # pixel size in arcsec (area per pixel = deltaPix**2)
        fwhm = 0.3  # full width half max of PSF

        # PSF specification

        kwargs_data = sim_util.data_configure_simple(numPix, deltaPix, exp_time, sigma_bkg)
        data_class = ImageData(**kwargs_data)
        sigma = util.fwhm2sigma(fwhm)
        x_grid, y_grid = util.make_grid(numPix=31, deltapix=0.05)
        from lenstronomy.LightModel.Profiles.gaussian import Gaussian
        gaussian = Gaussian()
        kernel_point_source = gaussian.function(x_grid, y_grid, amp=1., sigma=sigma, center_x=0, center_y=0)
        kernel_point_source /= np.sum(kernel_point_source)
        kernel_point_source = util.array2image(kernel_point_source)
        psf_error_map = np.zeros_like(kernel_point_source)
        self.kwargs_psf = {'psf_type': 'PIXEL', 'kernel_point_source': kernel_point_source,
                           'psf_error_map': psf_error_map}

        psf_class = PSF(**self.kwargs_psf)

        # '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)

if instrument_name == 'LSST':
        kwargs_instrument = LSST_camera
    else:
        raise ValueError("instrument name %s not supported! Chose among  %s " % (instrument_name, instrument_name_list))

    if observation_name == 'LSST_g_band':
        kwargs_observation = LSST_g_band_obs
    elif observation_name == 'LSST_r_band':
        kwargs_observation = LSST_r_band_obs
    elif observation_name == 'LSST_i_band':
        kwargs_observation = LSST_i_band_obs
    else:
        raise ValueError('observation name %s not supported! Chose among %s' %
                         (observation_name, observation_name_list))
    kwargs_data = util.merge_dicts(kwargs_instrument, kwargs_observation)
    return kwargs_data

:param x: x-coordinate in image plane
        :param y: y-coordinate in image plane
        :param theta_E: Einstein radius
        :param e1: eccentricity component
        :param e2: eccentricity component
        :param t: power law slope
        :param center_x: profile center
        :param center_y: profile center
        :return: f_xx, f_yy, f_xy
        """
        b, t, q, phi_G = self.param_conv(theta_E, e1, e2, t)
        # shift
        x_ = x - center_x
        y_ = y - center_y
        # rotate
        x__, y__ = util.rotate(x_, y_, phi_G)
        # evaluate
        f__xx, f__yy, f__xy = 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)
        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__
        f_xx = kappa + gamma1
        f_yy = kappa - gamma1
        f_xy = gamma2
        return f_xx, f_yy, f_xy

:param x: x-coordinate in image plane
        :param y: y-coordinate in image plane
        :param theta_E: Einstein radius
        :param e1: eccentricity component
        :param e2: eccentricity component
        :param t: power law slope
        :param center_x: profile center
        :param center_y: profile center
        :return: lensing potential
        """
        b, t, q, phi_G = self.param_conv(theta_E, e1, e2, t)
        # shift
        x_ = x - center_x
        y_ = y - center_y
        # rotate
        x__, y__ = util.rotate(x_, y_, phi_G)
        # evaluate
        f_ = self.epl_major_axis.function(x__, y__, b, t, q)
        # rotate back
        return f_

def source(self, numPix, deltaPix, center=None, image_orientation=True):
        """

        :param numPix: number of pixels per axes
        :param deltaPix: pixel size
        :param image_orientation: bool, if True, uses frame in orientation of the image, otherwise in RA-DEC coordinates
        :return: 2d surface brightness grid of the reconstructed source and Coordinates() instance of source grid
        """
        if image_orientation is True:
            Mpix2coord = self._coords.transform_pix2angle * deltaPix / self._deltaPix
            x_grid_source, y_grid_source = util.make_grid_transformed(numPix, Mpix2Angle=Mpix2coord)
            ra_at_xy_0, dec_at_xy_0 = x_grid_source[0], y_grid_source[0]
        else:
            x_grid_source, y_grid_source, ra_at_xy_0, dec_at_xy_0, x_at_radec_0, y_at_radec_0, Mpix2coord, Mcoord2pix = util.make_grid_with_coordtransform(
            numPix, deltaPix)

        center_x = 0
        center_y = 0
        if center is not None:
            center_x, center_y = center[0], center[1]
        elif len(self._kwargs_source_partial) > 0:
            center_x = self._kwargs_source_partial[0]['center_x']
            center_y = self._kwargs_source_partial[0]['center_y']
        x_grid_source += center_x
        y_grid_source += center_y

        coords_source = Coordinates(transform_pix2angle=Mpix2coord,

:param ny: number of pixels in y-axis
        :param transform_pix2angle: 2x2 matrix, mapping of pixel to coordinate
        :param ra_at_xy_0: ra coordinate at pixel (0,0)
        :param dec_at_xy_0: dec coordinate at pixel (0,0)
        :param supersampling_indexes: bool array of shape nx x ny, corresponding to pixels being super_sampled
        :param supersampling_factor: int, factor (per axis) of super-sampling
        :param flux_evaluate_indexes: bool array of shape nx x ny, corresponding to pixels being evaluated
        (for both low and high res). Default is None, replaced by setting all pixels to being evaluated.
        """
        super(AdaptiveGrid, self).__init__(transform_pix2angle, ra_at_xy_0, dec_at_xy_0)
        self._nx = nx
        self._ny = ny
        self._x_grid, self._y_grid = self.coordinate_grid(nx, ny)
        if flux_evaluate_indexes is None:
            flux_evaluate_indexes = np.ones_like(self._x_grid, dtype=bool)
        supersampled_indexes1d = util.image2array(supersampling_indexes)
        self._high_res_indexes1d = (supersampled_indexes1d) & (flux_evaluate_indexes)
        self._low_res_indexes1d = (np.invert(supersampled_indexes1d)) & (flux_evaluate_indexes)
        self._supersampling_factor = supersampling_factor
        self._num_sub = supersampling_factor * supersampling_factor
        self._x_low_res = self._x_grid[self._low_res_indexes1d]
        self._y_low_res = self._y_grid[self._low_res_indexes1d]
        self._num_low_res = len(self._x_low_res)