Not able to add NFW_MC lens to a mock image

[PrayaagKatta]

Hey, I was trying to create a mock image with the ImageModel class. It does work with a Pseudo Jaffe profile, but I am getting this error when I try with a NFW_MC profile -
`
File ~/.local/lib/python3.10/site-packages/lenstronomy/ImSim/image_model.py:333, in ImageModel.image(self, kwargs_lens, kwargs_source, kwargs_lens_light, kwargs_ps, kwargs_extinction, kwargs_special, unconvolved, source_add, lens_light_add, point_source_add)
316 def image(self, kwargs_lens=None, kwargs_source=None, kwargs_lens_light=None, kwargs_ps=None,
317 kwargs_extinction=None, kwargs_special=None, unconvolved=False, source_add=True, lens_light_add=True,
318 point_source_add=True):
319 """
320
321 make an image with a realisation of linear parameter values "param"
(...)
331 :return: 2d array of surface brightness pixels of the simulation
332 """
--> 333 return self._image(kwargs_lens, kwargs_source, kwargs_lens_light, kwargs_ps, kwargs_extinction, kwargs_special,
334 unconvolved, source_add, lens_light_add, point_source_add)

File ~/.local/lib/python3.10/site-packages/lenstronomy/ImSim/image_model.py:355, in ImageModel._image(self, kwargs_lens, kwargs_source, kwargs_lens_light, kwargs_ps, kwargs_extinction, kwargs_special, unconvolved, source_add, lens_light_add, point_source_add)
353 model = np.zeros(self.Data.num_pixel_axes)
354 if source_add is True:
--> 355 model += self._source_surface_brightness(kwargs_source, kwargs_lens, kwargs_extinction=kwargs_extinction,
356 kwargs_special=kwargs_special, unconvolved=unconvolved)
357 if lens_light_add is True:
358 model += self._lens_surface_brightness(kwargs_lens_light, unconvolved=unconvolved)

File ~/.local/lib/python3.10/site-packages/lenstronomy/ImSim/image_model.py:157, in ImageModel._source_surface_brightness(self, kwargs_source, kwargs_lens, kwargs_extinction, kwargs_special, unconvolved, de_lensed, k, update_pixelbased_mapping)
151 return self._source_surface_brightness_pixelbased(kwargs_source, kwargs_lens=kwargs_lens,
152 kwargs_extinction=kwargs_extinction,
153 kwargs_special=kwargs_special,
154 unconvolved=unconvolved, de_lensed=de_lensed, k=k,
155 update_mapping=update_pixelbased_mapping)
156 else:
--> 157 return self._source_surface_brightness_analytical(kwargs_source, kwargs_lens=kwargs_lens,
158 kwargs_extinction=kwargs_extinction,
159 kwargs_special=kwargs_special,
160 unconvolved=unconvolved, de_lensed=de_lensed, k=k)

File ~/.local/lib/python3.10/site-packages/lenstronomy/ImSim/image_model.py:180, in ImageModel._source_surface_brightness_analytical(self, kwargs_source, kwargs_lens, kwargs_extinction, kwargs_special, unconvolved, de_lensed, k)
178 source_light = self.SourceModel.surface_brightness(ra_grid, dec_grid, kwargs_source, k=k)
179 else:
--> 180 source_light = self.source_mapping.image_flux_joint(ra_grid, dec_grid, kwargs_lens, kwargs_source, k=k)
181 source_light *= self._extinction.extinction(ra_grid, dec_grid, kwargs_extinction=kwargs_extinction,
182 kwargs_special=kwargs_special)
184 # multiply with primary beam before convolution

File ~/.local/lib/python3.10/site-packages/lenstronomy/ImSim/image2source_mapping.py:120, in Image2SourceMapping.image_flux_joint(self, x, y, kwargs_lens, kwargs_source, k)
110 """
111
112 :param x: coordinate in image plane
(...)
117 :return: surface brightness of all joint light components at image position (x, y)
118 """
119 if self._multi_source_plane is False:
--> 120 x_source, y_source = self._lensModel.ray_shooting(x, y, kwargs_lens)
121 return self._lightModel.surface_brightness(x_source, y_source, kwargs_source, k=k)
122 else:

File ~/.local/lib/python3.10/site-packages/lenstronomy/LensModel/lens_model.py:117, in LensModel.ray_shooting(self, x, y, kwargs, k)
105 def ray_shooting(self, x, y, kwargs, k=None):
106 """
107 maps image to source position (inverse deflection)
108
(...)
115 :return: source plane positions corresponding to (x, y) in the image plane
116 """
--> 117 return self.lens_model.ray_shooting(x, y, kwargs, k=k)

File ~/.local/lib/python3.10/site-packages/lenstronomy/LensModel/single_plane.py:26, in SinglePlane.ray_shooting(self, x, y, kwargs, k)
14 def ray_shooting(self, x, y, kwargs, k=None):
15 """
16 maps image to source position (inverse deflection)
17 :param x: x-position (preferentially arcsec)
(...)
23 :return: source plane positions corresponding to (x, y) in the image plane
24 """
---> 26 dx, dy = self.alpha(x, y, kwargs, k=k)
27 return x - dx, y - dy

File ~/.local/lib/python3.10/site-packages/lenstronomy/LensModel/single_plane.py:92, in SinglePlane.alpha(self, x, y, kwargs, k)
90 for i, func in enumerate(self.func_list):
91 if bool_list[i] is True:
---> 92 f_x_i, f_y_i = func.derivatives(x, y, **kwargs[i])
93 f_x += f_x_i
94 f_y += f_y_i

File ~/.local/lib/python3.10/site-packages/lenstronomy/LensModel/Profiles/nfw_mass_concentration.py:99, in NFWMC.derivatives(self, x, y, logM, concentration, center_x, center_y)
95 def derivatives(self, x, y, logM, concentration, center_x=0, center_y=0):
96 """
97 returns df/dx and df/dy of the function (integral of NFW)
98 """
---> 99 Rs, alpha_Rs = self._m_c2deflections(logM, concentration)
100 return self._nfw.derivatives(x, y, Rs, alpha_Rs, center_x, center_y)

File ~/.local/lib/python3.10/site-packages/lenstronomy/LensModel/Profiles/nfw_mass_concentration.py:54, in NFWMC._m_c2deflections(self, logM, concentration)
52 return self._Rs_static, self._alpha_Rs_static
53 M = 10 ** logM
---> 54 Rs, alpha_Rs = self._lens_cosmo.nfw_physical2angle(M, concentration)
55 return Rs, alpha_Rs

File ~/.local/lib/python3.10/site-packages/lenstronomy/Cosmo/lens_cosmo.py:195, in LensCosmo.nfw_physical2angle(self, M, c)
187 def nfw_physical2angle(self, M, c):
188 """
189 converts the physical mass and concentration parameter of an NFW profile into the lensing quantities
190
(...)
193 :return: Rs_angle (angle at scale radius) (in units of arcsec), alpha_Rs (observed bending angle at the scale radius
194 """
--> 195 rho0, Rs, r200 = self.nfwParam_physical(M, c)
196 Rs_angle = Rs / self.dd / const.arcsec # Rs in arcsec
197 alpha_Rs = rho0 * (4 * Rs ** 2 * (1 + np.log(1. / 2.)))

File ~/.local/lib/python3.10/site-packages/lenstronomy/Cosmo/lens_cosmo.py:208, in LensCosmo.nfwParam_physical(self, M, c)
200 def nfwParam_physical(self, M, c):
201 """
202 returns the NFW parameters in physical units
203
(...)
206 :return: rho0 [Msun/Mpc^3], Rs [Mpc], r200 [Mpc]
207 """
--> 208 r200 = self.nfw_param.r200_M(M * self.h, self.z_lens) / self.h # physical radius r200
209 rho0 = self.nfw_param.rho0_c(c, self.z_lens) * self.h2 # physical density in M_sun/Mpc3
210 Rs = r200/c

File ~/.local/lib/python3.10/site-packages/lenstronomy/Cosmo/nfw_param.py:64, in NFWParam.r200_M(self, M, z)
54 def r200_M(self, M, z):
55 """
56 computes the radius R_200 crit of a halo of mass M in physical mass M/h
57
(...)
62 :return: radius R_200 in physical Mpc/h
63 """
---> 64 return (3M/(4np.pi*self.rhoc_z(z)*200))**(1./3.)

File ~/.local/lib/python3.10/site-packages/lenstronomy/Cosmo/nfw_param.py:33, in NFWParam.rhoc_z(self, z)
27 def rhoc_z(self, z):
28 """
29
30 :param z: redshift
31 :return: critical density of the universe at redshift z in physical units [h^2 M_sun Mpc^-3]
32 """
---> 33 return self.rhoc * (self.cosmo.efunc(z)) ** 2

File ~/.local/lib/python3.10/site-packages/astropy/cosmology/flrw/lambdacdm.py:734, in FlatLambdaCDM.efunc(self, z)
727 # We override this because it takes a particularly simple
728 # form for a cosmological constant
729 Or = self._Ogamma0 + (
730 self._Onu0
731 if not self._massivenu
732 else self._Ogamma0 * self.nu_relative_density(z)
733 )
--> 734 zp1 = aszarr(z) + 1.0 # (converts z [unit] -> z [dimensionless])
736 return np.sqrt(zp1**3 * (Or * zp1 + self._Om0) + self._Ode0)

File ~/.local/lib/python3.10/site-packages/astropy/cosmology/utils.py:149, in aszarr(z)
147 return z
148 # not one of the preferred types: Number / array ducktype
--> 149 return Quantity(z, cu.redshift).value

File ~/.local/lib/python3.10/site-packages/astropy/units/quantity.py:552, in Quantity.new(cls, value, unit, dtype, copy, order, subok, ndmin)
548 # check that array contains numbers or long int objects
549 if value.dtype.kind in "OSU" and not (
550 value.dtype.kind == "O" and isinstance(value.item(0), numbers.Number)
551 ):
--> 552 raise TypeError("The value must be a valid Python or Numpy numeric type.")
554 # by default, cast any integer, boolean, etc., to float
555 if float_default and value.dtype.kind in "iuO":

TypeError: The value must be a valid Python or Numpy numeric type.
`
Thanks in advance for your help!

[sibirrer]

@PrayaagKatta is this still relevant? Just checking. What astropy version are you using?
Simon

[PrayaagKatta]

Hey @sibirrer, yes I am still facing this issue. If it helps, I can also give the function I am using to make the mock, let me know! The astropy version is 5.3.1

[sibirrer]

is the input mass a float and the redshift as well? Please state the entire function and input that leads to this error message

[PrayaagKatta]

The function is -

def add_perturber_to_mock(perturber_log_mass,perturber_redshift,perturber_type,perturber_coords_list,
                          Ra, Rs, perturber_conc, source_redshift):
    #Defining everything needed to make mock image first
    # import the PixelGrid() class #
    deltaPix = 0.01  # size of pixel in angular coordinates #
    numPix = 400
    # setup the keyword arguments to create the Data() class #
    ra_at_xy_0, dec_at_xy_0 = -numPix*deltaPix/2.0, -numPix*deltaPix/2.0 # coordinate in angles (RA/DEC) at the position of the pixel edge (0,0)

    # generate the coordinate grid and image properties (we only read out the relevant lines we need)
    _, _, ra_at_xy_0, dec_at_xy_0, _, _, Mpix2coord, _ = util.make_grid_with_coordtransform(numPix=numPix, deltapix=deltaPix, center_ra=0, center_dec=0, subgrid_res=1, inverse=False)
    coords = Coordinates(Mpix2coord, ra_at_xy_0, dec_at_xy_0)

    # transform_pix2angle = np.array([[1, 0], [0, 1]]) * deltaPix  # linear translation matrix of a shift in pixel in a shift in coordinates
    kwargs_pixel = {'nx': numPix, 'ny': numPix,  # 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': Mpix2coord} 

    pixel_grid = PixelGrid(**kwargs_pixel)
    # return the list of pixel coordinates #
    x_coords, y_coords = pixel_grid.pixel_coordinates

    # import the PSF() class #
    from lenstronomy.Data.psf import PSF
    kwargs_psf = {'psf_type': 'GAUSSIAN',  # type of PSF model (supports 'GAUSSIAN' and 'PIXEL')
                'fwhm': 0.05,  # full width at half maximum of the Gaussian PSF (in angular units)
                'pixel_size': deltaPix  # angular scale of a pixel (required for a Gaussian PSF to translate the FWHM into a pixel scale)
                }
    psf = PSF(**kwargs_psf)
    # return the pixel kernel corresponding to a point source # 
    kernel = psf.kernel_point_source

    # define the numerics 
    kwargs_numerics = {'supersampling_factor': 3, 'supersampling_convolution': False}


    previous_kwargs_result = {'kwargs_lens': [{'theta_E': 1.180804036937751, 'gamma': 2.0, 'e1': 0.2543259566546414, 'e2': 0.4228806445841423, 'center_x': 0.13179582388628278, 'center_y': -0.16449843135712475}, 
                                              {'gamma1': 0.08771548427110364, 'gamma2': 0.043840524845854156, 'ra_0': 0.11308158275357881, 'dec_0': -0.14766238077236038},
                                             ], 
                            'kwargs_source': [{'amp': 811.5505715444301, 'R_sersic': 0.27090161662895057, 'n_sersic': 3.357015167672101, 'e1': -0.13767048426057404, 'e2': -0.3228695558557331, 'center_x': -0.09685991514216577, 'center_y': 0.17716316728288928}], 
                            'kwargs_lens_light': [{'amp': 1756.6771358565247, 'R_sersic': 0.645530455839, 'n_sersic': 0.3210387621003749, 'e1': 0.042091786561357286, 'e2': -0.05018194983550343, 'center_x': 0.264383495143268, 'center_y': -0.2540911599976115}, 
                                                  {'amp': 7795.693094849382, 'R_sersic': 0.16743102249165667, 'n_sersic': 1.1223134912522115, 'e1': 0.0556364725508792, 'e2': -0.0024589201740620875, 'center_x': 0.004994380615012455, 'center_y': 0.007059840518109048}], 
                            'kwargs_ps': [], 'kwargs_special': {}, 'kwargs_extinction': []}

    lens_redshift_list=[0.5,0.5]
    lens_model_list = ['EPL', 'SHEAR']
    source_model_list = ['SERSIC_ELLIPSE']
    lens_light_model_list = ['SERSIC_ELLIPSE','SERSIC_ELLIPSE']

    ######################################################################
    #Including perturber stuff
    if(perturber_type == "PJAFFE"):
        cosmo = default_cosmology.get()
        lensCosmo = LensCosmo(z_lens=perturber_redshift, z_source=source_redshift, cosmo=cosmo)
        Sigma_crit = lensCosmo.sigma_crit_angle
        Sigma0 = (10**perturber_log_mass)/(2*np.pi*Ra*Rs)
        sigma0 = Sigma0/Sigma_crit
        previous_kwargs_result['kwargs_lens'].append({'sigma0': sigma0, 'Ra': Ra, 'Rs': Rs, 'center_x': perturber_coords[0], 'center_y': perturber_coords[1]})
    elif(perturber_type == "NFW_MC"):
        previous_kwargs_result['kwargs_lens'].append({'logM': perturber_log_mass, 'concentration': perturber_conc, 'center_x': perturber_coords[0], 'center_y': perturber_coords[1]})

    lens_redshift_list.append(perturber_redshift)
    lens_model_list.append(perturber_type)
    ######################################################################
    lensModel = LensModel(lens_model_list=lens_model_list)
    lightModel_source = LightModel(light_model_list=source_model_list)
    lightModel_lens = LightModel(light_model_list=lens_light_model_list)

    # initialize the Image model class by combining the modules we created above #
    imageModel = ImageModel(data_class=pixel_grid, psf_class=psf, lens_model_class=lensModel,
                            source_model_class=lightModel_source,
                            lens_light_model_class=lightModel_lens,
                            point_source_class=None, # in this example, we do not simulate point source.
                            kwargs_numerics=kwargs_numerics)
    # simulate image with the parameters we have defined above #
    image = imageModel.image(kwargs_lens=previous_kwargs_result['kwargs_lens'], kwargs_source=previous_kwargs_result['kwargs_source'],
                            kwargs_lens_light=previous_kwargs_result['kwargs_lens_light'], kwargs_ps=previous_kwargs_result['kwargs_ps'])
    
    # # Adding noise
    exp_time = 180*38 
    background_rms = 0.0615377  # background rms value from data
    poisson = image_util.add_poisson(image, exp_time=exp_time)
    bkg = image_util.add_background(image, sigma_bkd=background_rms)
    image_noisy = image + bkg + poisson

    ##########################################################################
    lightModel_lens_without_lenslight = LightModel(light_model_list=[])

    # initialize the Image model class by combining the modules we created above #
    imageModel_without_lenslight = ImageModel(data_class=pixel_grid, psf_class=psf, lens_model_class=lensModel,
                                              source_model_class=lightModel_source,
                                              lens_light_model_class=lightModel_lens_without_lenslight,
                                              point_source_class=None, # in this example, we do not simulate point source.
                                              kwargs_numerics=kwargs_numerics)
    # simulate image with the parameters we have defined above #
    image_without_lenslight = imageModel_without_lenslight.image(kwargs_lens=previous_kwargs_result['kwargs_lens'], kwargs_source=previous_kwargs_result['kwargs_source'],
                                                                 kwargs_lens_light=[{}], kwargs_ps=previous_kwargs_result['kwargs_ps'])
    
    # Adding noise
    image_noisy_without_lenslight = image_without_lenslight + bkg + poisson
    f, axes = plt.subplots(1, 3, figsize=(16, 8), sharex=False, sharey=False)
    axes[0].matshow(np.log10(image), origin='lower')
    axes[0].set_title("Model")
    axes[1].matshow(np.log10(image_noisy), origin='lower')
    axes[1].set_title("Model+Noise")
    axes[2].matshow(np.log10(image_noisy_without_lenslight), origin='lower')
    axes[2].set_title("Model without lens light")
    f.tight_layout()
    f.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.1, hspace=None)
    if(perturber_type=='PJAFFE'):
        f.suptitle('PJ_mass-{}_Ra-{}_Rs-{}'.format(perturber_log_mass,Ra,Rs))
    elif(perturber_type=='NFW_MC'):
        f.suptitle('NFW_mass-{}_concentration_{}'.format(perturber_log_mass,perturber_conc))
    plt.show()
    


    kwargs_data = {'image_data': image_noisy_without_lenslight,
                   'noise_map': bkg + poisson,
                   'ra_at_xy_0': ra_at_xy_0,
                   'dec_at_xy_0': dec_at_xy_0,
                   'transform_pix2angle': Mpix2coord
                  }

    return image_noisy, image_noisy_without_lenslight, kwargs_data

And the input is -

#Setting up perturbers

perturber_redshift = 0.5
source_redshift = 2
perturber_log_mass = 10
perturber_type = "NFW_MC"  #"NFW_MC" or "PJAFFE"
perturber_coords = [perturber_coords_list[0][0], perturber_coords_list[1][0]]
Ra = 0.001
Rs = 0.166
perturber_conc = 2

mock_image_with_lens_light, mock_image, kwargs_data = add_perturber_to_mock(perturber_log_mass,perturber_redshift,perturber_type,perturber_coords, 
                                                                            Ra, Rs, perturber_conc, source_redshift)

[sibirrer]

hmm, unclear what's happening here. Perhaps @dangilman knows?