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

def test_delete_interpol_caches(self):
        x, y = util.make_grid(numPix=20, deltapix=1.)
        gauss = Gaussian()
        flux = gauss.function(x, y, amp=1., center_x=0., center_y=0., sigma=1.)
        image = util.array2image(flux)

        light_model_list = ['INTERPOL', 'INTERPOL']
        kwargs_list = [
            {'image': image, 'scale': 1, 'phi_G': 0, 'center_x': 0, 'center_y': 0},
            {'image': image, 'scale': 1, 'phi_G': 0, 'center_x': 0, 'center_y': 0}
        ]
        lightModel = LightModel(light_model_list=light_model_list)
        output = lightModel.surface_brightness(x, y, kwargs_list)
        for func in lightModel.func_list:
            assert hasattr(func, '_image_interp')
        lightModel.delete_interpol_caches()
        for func in lightModel.func_list:
            assert not hasattr(func, '_image_interp')

def setup(self):
        self.supersampling_factor = 3
        lightModel = LightModel(light_model_list=['GAUSSIAN'])
        self.delta_pix = 1.
        x, y = util.make_grid(20, deltapix=self.delta_pix)
        x_sub, y_sub = util.make_grid(20*self.supersampling_factor, deltapix=self.delta_pix/self.supersampling_factor)
        kwargs = [{'amp': 1, 'sigma': 2, 'center_x': 0, 'center_y': 0}]
        flux = lightModel.surface_brightness(x, y, kwargs)
        self.model = util.array2image(flux)
        flux_sub = lightModel.surface_brightness(x_sub, y_sub, kwargs)
        self.model_sub = util.array2image(flux_sub)

        x, y = util.make_grid(5, deltapix=self.delta_pix)
        kwargs_kernel = [{'amp': 1, 'sigma': 1, 'center_x': 0, 'center_y': 0}]
        kernel = lightModel.surface_brightness(x, y, kwargs_kernel)
        self.kernel = util.array2image(kernel) / np.sum(kernel)

        x_sub, y_sub = util.make_grid(5*self.supersampling_factor, deltapix=self.delta_pix/self.supersampling_factor)
        kernel_sub = lightModel.surface_brightness(x_sub, y_sub, kwargs_kernel)
        self.kernel_sub = util.array2image(kernel_sub) / np.sum(kernel_sub)

x, y = util.make_grid(20, deltapix=self.delta_pix)
        x_sub, y_sub = util.make_grid(20*self.supersampling_factor, deltapix=self.delta_pix/self.supersampling_factor)
        kwargs = [{'amp': 1, 'sigma': 2, 'center_x': 0, 'center_y': 0}]
        flux = lightModel.surface_brightness(x, y, kwargs)
        self.model = util.array2image(flux)
        flux_sub = lightModel.surface_brightness(x_sub, y_sub, kwargs)
        self.model_sub = util.array2image(flux_sub)

        x, y = util.make_grid(5, deltapix=self.delta_pix)
        kwargs_kernel = [{'amp': 1, 'sigma': 1, 'center_x': 0, 'center_y': 0}]
        kernel = lightModel.surface_brightness(x, y, kwargs_kernel)
        self.kernel = util.array2image(kernel) / np.sum(kernel)

        x_sub, y_sub = util.make_grid(5*self.supersampling_factor, deltapix=self.delta_pix/self.supersampling_factor)
        kernel_sub = lightModel.surface_brightness(x_sub, y_sub, kwargs_kernel)
        self.kernel_sub = util.array2image(kernel_sub) / np.sum(kernel_sub)

def test_interp_func_scaled(self):

        numPix = 101
        deltaPix = 0.1
        x_grid_interp, y_grid_interp = util.make_grid(numPix,deltaPix)
        sis = SIS()
        kwargs_SIS = {'theta_E':1., 'center_x': 0.5, 'center_y': -0.5}
        f_sis = sis.function(x_grid_interp, y_grid_interp, **kwargs_SIS)
        f_x_sis, f_y_sis = sis.derivatives(x_grid_interp, y_grid_interp, **kwargs_SIS)
        f_xx_sis, f_yy_sis, f_xy_sis = sis.hessian(x_grid_interp, y_grid_interp, **kwargs_SIS)
        x_axes, y_axes = util.get_axes(x_grid_interp, y_grid_interp)
        kwargs_interp = {'grid_interp_x': x_axes, 'grid_interp_y': y_axes, 'f_': util.array2image(f_sis), 'f_x': util.array2image(f_x_sis), 'f_y': util.array2image(f_y_sis), 'f_xx': util.array2image(f_xx_sis), 'f_yy': util.array2image(f_yy_sis), 'f_xy': util.array2image(f_xy_sis)}
        interp_func = InterpolScaled(grid=False)
        x, y = 1., 1.
        alpha_x, alpha_y = interp_func.derivatives(x, y, scale_factor=1, **kwargs_interp)
        assert alpha_x == 0.31622776601683794

        f_ = interp_func.function(x, y, scale_factor=1., **kwargs_interp)
        npt.assert_almost_equal(f_, 1.5811388300841898)

        f_xx, f_yy, f_xy = interp_func.hessian(x, y, scale_factor=1., **kwargs_interp)
        npt.assert_almost_equal(f_xx, 0.56920997883030822, decimal=8)
        npt.assert_almost_equal(f_yy, 0.063245553203367583, decimal=8)
        npt.assert_almost_equal(f_xy, -0.18973665961010275, decimal=8)

        x_grid, y_grid = util.make_grid(10, deltaPix)
        f_xx, f_yy, f_xy = interp_func.hessian(x_grid, y_grid, scale_factor=1., **kwargs_interp)
        npt.assert_almost_equal(f_xx[0], 0, decimal=2)

def pixel_kernel(self, num_pix):
        """
        computes a pixelized kernel from the MGE parameters

        :param num_pix: int, size of kernel (odd number per axis)
        :return: pixel kernel centered
        """
        from lenstronomy.LightModel.Profiles.gaussian import MultiGaussian
        mg = MultiGaussian()
        x, y = util.make_grid(numPix=num_pix, deltapix=self._pixel_scale)
        kernel = mg.function(x, y, amp=self._fraction_list, sigma=self._sigmas_scaled)
        kernel = util.array2image(kernel)
        return kernel / np.sum(kernel)

:param psf_type:
        :param fwhm:
        :param pixel_grid:
        :return:
        """
        # psf_type: 'NONE', 'gaussian', 'pixel'
        # 'pixel': kernel, kernel_large
        # 'gaussian': 'sigma', 'truncate'
        if psf_type == 'GAUSSIAN':
            sigma = util.fwhm2sigma(fwhm)
            sigma_axis = sigma
            x_grid, y_grid = util.make_grid(kernelsize, deltaPix)
            kernel_large = self.gaussian.function(x_grid, y_grid, amp=1., sigma_x=sigma_axis, sigma_y=sigma_axis, center_x=0, center_y=0)
            kernel_large /= np.sum(kernel_large)
            kernel_large = util.array2image(kernel_large)
            kernel_pixel = kernel_util.pixel_kernel(kernel_large)
            kwargs_psf = {'psf_type': psf_type, 'fwhm': fwhm, 'truncation': truncate*fwhm, 'kernel_point_source': kernel_large, 'kernel_pixel': kernel_pixel, 'pixel_size': deltaPix}
        elif psf_type == 'PIXEL':
            kernel_large = copy.deepcopy(kernel)
            kernel_large = kernel_util.cut_psf(kernel_large, psf_size=kernelsize)
            kernel_small = copy.deepcopy(kernel)
            kernel_small = kernel_util.cut_psf(kernel_small, psf_size=kernelsize)
            #kwargs_psf = {'psf_type': "PIXEL", 'kernel_pixel': kernel_small, 'kernel_point_source': kernel_large}
            kwargs_psf = {'psf_type': "PIXEL", 'kernel_point_source': kernel_large}
        elif psf_type == 'NONE':
            kwargs_psf = {'psf_type': 'NONE'}
        else:
            raise ValueError("psf type %s not supported!" % psf_type)
        return kwargs_psf

x_axes_sub, y_axes_sub = util.get_axes(x_grid_sub, y_grid_sub)
    from lenstronomy.LensModel.Profiles.interpol import Interpol
    interp_func = Interpol()
    interp_func.do_interp(x_axes_sub, y_axes_sub, f_sub, f_x_sub, f_y_sub)
    # compute lensing quantities on sparser grid
    x_axes, y_axes = util.get_axes(x_grid, y_grid)
    f_ = interp_func.function(x_grid, y_grid)
    f_x, f_y = interp_func.derivatives(x_grid, y_grid)
    # numerical differentials for second order differentials
    from lenstronomy.LensModel.lens_model import LensModel
    lens_model = LensModel(lens_model_list=['INTERPOL'])
    kwargs = [{'grid_interp_x': x_axes_sub, 'grid_interp_y': y_axes_sub, 'f_': f_sub,
               'f_x': f_x_sub, 'f_y': f_y_sub}]
    f_xx, f_xy, f_yx, f_yy = lens_model.hessian(x_grid, y_grid, kwargs, diff=0.00001)
    kwargs_interpol = {'grid_interp_x': x_axes, 'grid_interp_y': y_axes, 'f_': util.array2image(f_),
                       'f_x': util.array2image(f_x), 'f_y': util.array2image(f_y), 'f_xx': util.array2image(f_xx),
                       'f_xy': util.array2image(f_xy), 'f_yy': util.array2image(f_yy)}
    return kwargs_interpol

:return:
        """
        x_s, y_s = self.mapping_IS(x_grid, y_grid, kwargs_else, **kwargs_lens)


        image = np.zeros(numPix**2)
        center_x = kwargs_source['center_x']
        center_y = kwargs_source['center_y']
        num_param_shapelets = (num_order + 2) * (num_order + 1) / 2
        H_x, H_y = self.shapelets.pre_calc(x_s, y_s, beta, num_order, center_x, center_y)
        n1 = 0
        n2 = 0
        for i in range(len(param) - num_param_shapelets, len(param)):
            kwargs_source_shapelet = {'center_x': center_x, 'center_y': center_y, 'n1': n1, 'n2': n2, 'beta': beta, 'amp': param[i]}
            image_i = self.shapelets.function(H_x, H_y, **kwargs_source_shapelet)
            image_i = util.array2image(image_i)
            image_i = self.re_size_convolve(image_i, numPix, deltaPix, subgrid_res, kwargs_psf)
            image += util.image2array(image_i*mask)
            if n1 == 0:
                n1 = n2 + 1
                n2 = 0
            else:
                n1 -= 1
                n2 += 1
        return image

y_shear = y_grid - f_y_shear

    if foreground_shear:
        f_x_shear1, f_y_shear1 = shear.derivatives(x_grid, y_grid, e1=kwargs_else['gamma1_foreground']*strength_multiply, e2=kwargs_else['gamma2_foreground']*strength_multiply)
    else:
        f_x_shear1, f_y_shear1 = 0, 0
    x_foreground = x_grid - f_x_shear1
    y_foreground = y_grid - f_y_shear1

    center_x = np.mean(x_grid)
    center_y = np.mean(y_grid)
    radius = (np.max(x_grid) - np.min(x_grid))/4
    circle_shear = util_maskl.mask_sphere(x_shear, y_shear, center_x, center_y, radius)
    circle_foreground = util_maskl.mask_sphere(x_foreground, y_foreground, center_x, center_y, radius)
    f, ax = plt.subplots(1, 1, figsize=(16, 8), sharex=False, sharey=False)
    im = ax.matshow(util.array2image(np.log10(kwargs_data['image_data'])), origin='lower', alpha=0.5)
    im = ax.matshow(util.array2image(circle_shear), origin='lower', alpha=0.5, cmap="jet")
    im = ax.matshow(util.array2image(circle_foreground), origin='lower', alpha=0.5)
    #f.show()
    return f, ax

:param smoothing_scale: float or None, Gaussian FWHM of a smoothing kernel applied before plotting
    :return: matplotlib instance with different panels
    """
    kwargs_grid = sim_util.data_configure_simple(num_pix, delta_pix, center_ra=center_ra, center_dec=center_dec)
    _coords = ImageData(**kwargs_grid)
    _frame_size = num_pix * delta_pix
    ra_grid, dec_grid = _coords.pixel_coordinates

    extensions = LensModelExtensions(lensModel=lensModel)
    ra_grid1d = util.image2array(ra_grid)
    dec_grid1d = util.image2array(dec_grid)
    lambda_rad, lambda_tan, orientation_angle, dlambda_tan_dtan, dlambda_tan_drad, dlambda_rad_drad, dlambda_rad_dtan, dphi_tan_dtan, dphi_tan_drad, dphi_rad_drad, dphi_rad_dtan = extensions.radial_tangential_differentials(
        ra_grid1d, dec_grid1d, kwargs_lens=kwargs_lens, center_x=center_ra, center_y=center_dec, smoothing_3rd=differential_scale, smoothing_2nd=None)

    lambda_rad2d, lambda_tan2d, orientation_angle2d, dlambda_tan_dtan2d, dlambda_tan_drad2d, dlambda_rad_drad2d, dlambda_rad_dtan2d, dphi_tan_dtan2d, dphi_tan_drad2d, dphi_rad_drad2d, dphi_rad_dtan2d = util.array2image(lambda_rad), \
                                            util.array2image(lambda_tan), util.array2image(orientation_angle), util.array2image(dlambda_tan_dtan), util.array2image(dlambda_tan_drad), util.array2image(dlambda_rad_drad), util.array2image(dlambda_rad_dtan), \
                                            util.array2image(dphi_tan_dtan), util.array2image(dphi_tan_drad), util.array2image(dphi_rad_drad), util.array2image(dphi_rad_dtan)

    if smoothing_scale is not None:
        lambda_rad2d = ndimage.gaussian_filter(lambda_rad2d, sigma=smoothing_scale/delta_pix)
        dlambda_rad_drad2d = ndimage.gaussian_filter(dlambda_rad_drad2d, sigma=smoothing_scale/delta_pix)
        lambda_tan2d = np.abs(lambda_tan2d)
        # the magnification cut is made to make a stable integral/convolution
        lambda_tan2d[lambda_tan2d > 100] = 100
        lambda_tan2d = ndimage.gaussian_filter(lambda_tan2d, sigma=smoothing_scale/delta_pix)
        # the magnification cut is made to make a stable integral/convolution
        dlambda_tan_dtan2d[dlambda_tan_dtan2d > 100] = 100
        dlambda_tan_dtan2d[dlambda_tan_dtan2d < -100] = -100
        dlambda_tan_dtan2d = ndimage.gaussian_filter(dlambda_tan_dtan2d, sigma=smoothing_scale/delta_pix)
        orientation_angle2d = ndimage.gaussian_filter(orientation_angle2d, sigma=smoothing_scale/delta_pix)
        dphi_tan_dtan2d = ndimage.gaussian_filter(dphi_tan_dtan2d, sigma=smoothing_scale/delta_pix)