Skip to content

likelihood

photXclus.likelihood

NFW_r(r, rx, integr_c)

Evaluate the Navarro-Frank-White (NFW) model at given physical radii

Parameters:

Name Type Description Default
r array_like

Array of physical radii

 required
rx float

Scale radius of the model

 required
integr_c float

Integral of the model

 required

Returns:

Name Type Description
NFW_model array_like

Array of the Navarro-Frank-White model evaluated at the input physical radii
Source code in photXclus/model.py 
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
@nb.njit()
def NFW_r(r, rx, integr_c):
    """
    Evaluate the Navarro-Frank-White (NFW) model at given physical radii

    Parameters
    ----------
    r : array_like
        Array of physical radii
    rx : float
        Scale radius of the model
    integr_c : float
        Integral of the model

    Returns
    -------
    NFW_model : array_like
        Array of the Navarro-Frank-White model evaluated at the input physical radii
    """
    rs = rx/3.
    return integr_c * psi_r_cut(r, rs)

clus_field_zr(z, r, ngal_c, muz_c, bin_zc, sigz_z, ngal_f, binz_ampl_f, bin_z, C_Omega, rx, deltac)

Evaluate the total model at given redshifts and physical radii

Parameters:

Name Type Description Default
z array_like

Array of redshifts

 required
r array_like

Array of physical radii

 required
ngal_c float

Number of galaxies in the cluster

 required
muz_c float

Mean of the cluster redshift distribution

 required
bin_zc array_like

Array of redshift bins

 required
sigz_z array_like

Array of standard deviation of the galaxy photometric redshift (sigma_zphot(z))

 required
ngal_f float

Number density of galaxies in the field in degrees^-2

 required
binz_ampl_f array_like

Amplitude of the field model at each redshift bin

 required
bin_z array_like

Array of redshift bins

 required
C_Omega float

Area of the field of view in degrees^-2

 required
rx float

Scale radius of the cluster model

 required
deltac float

Delta parameter for the field model

 required

Returns:

Name Type Description
total_model array_like

Array of the total model evaluated at the input redshifts and physical radii
Source code in photXclus/model.py 
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
@nb.njit()
def clus_field_zr(z, r, ngal_c, 
                  muz_c, bin_zc, sigz_z, 
                  ngal_f, 
                  binz_ampl_f, bin_z, C_Omega, rx, deltac):
    """
    Evaluate the total model at given redshifts and physical radii

    Parameters
    ----------
    z : array_like
        Array of redshifts
    r : array_like
        Array of physical radii
    ngal_c : float
        Number of galaxies in the cluster
    muz_c : float
        Mean of the cluster redshift distribution
    bin_zc : array_like
        Array of redshift bins
    sigz_z : array_like
        Array of standard deviation of the galaxy photometric redshift (sigma_zphot(z))
    ngal_f : float
        Number density of galaxies in the field in degrees^-2
    binz_ampl_f : array_like
        Amplitude of the field model at each redshift bin
    bin_z : array_like
        Array of redshift bins
    C_Omega : float
        Area of the field of view in degrees^-2
    rx : float
        Scale radius of the cluster model
    deltac : float
        Delta parameter for the field model

    Returns
    -------
    total_model : array_like
        Array of the total model evaluated at the input redshifts and physical radii
    """
    if deltac == 1:
        pclus = clus_z(z, muz_c, bin_zc, sigz_z) * clus_r(r, rx)
        pfield = field_z(z, binz_ampl_f, bin_z, sigz_z) * field_r(r, C_Omega)
        ptot = ngal_c * pclus + C_Omega * ngal_f * pfield
    else:
        pfield = field_z(z, binz_ampl_f, bin_z, sigz_z)
        ptot = C_Omega * ngal_f * pfield
    return ptot

clus_r(r, rx, integr_c=1.0, model='Plummer')

Evaluate the cluster radial model at given physical radii

Parameters:

Name Type Description Default
r array_like

Array of physical radii

 required
rx float

Scale radius of the cluster model

 required
integr_c float

Integral of the cluster model, used for the NFW model

 1.0
model str

Model to be used, either 'Plum' or 'NFW'

 'Plummer'

Returns:

Name Type Description
clus_model array_like

Array of the cluster model evaluated at the input physical radii
Source code in photXclus/model.py 
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
@nb.njit()
def clus_r(r, rx, integr_c = 1., model='Plummer'):
    """
    Evaluate the cluster radial model at given physical radii

    Parameters
    ----------
    r : array_like
        Array of physical radii
    rx : float
        Scale radius of the cluster model
    integr_c : float, optional
        Integral of the cluster model, used for the NFW model
    model : str, optional
        Model to be used, either 'Plum' or 'NFW'

    Returns
    -------
    clus_model : array_like
        Array of the cluster model evaluated at the input physical radii
    """
    if model == 'Plummer':
        return plum_r(r, rx)
    elif model == 'NFW':
        return NFW_r(r, rx, integr_c)
    else:
        raise NotImplementedError(f"{model} radial model not implemented")

clus_z(z, muz_c, z_binc, sigz_z)

Evaluate the cluster redshift model at given redshifts

Parameters:

Name Type Description Default
z array_like

Array of redshifts

 required
muz_c float

Mean of the cluster redshift distribution

 required
z_binc array_like

Array of redshift bins

 required
sigz_z array_like

Array of standard deviation of the galaxy photometric redshift (sigma_zphot(z))

 required

Returns:

Name Type Description
clus_model array_like

Array of the cluster model evaluated at the input redshifts
Source code in photXclus/model.py 
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
@nb.njit()
def clus_z(z, muz_c, z_binc, sigz_z ):
    """
    Evaluate the cluster redshift model at given redshifts

    Parameters
    ----------
    z : array_like
        Array of redshifts
    muz_c : float
        Mean of the cluster redshift distribution
    z_binc : array_like
        Array of redshift bins
    sigz_z : array_like
        Array of standard deviation of the galaxy photometric redshift (sigma_zphot(z))

    Returns
    -------
    clus_model : array_like
        Array of the cluster model evaluated at the input redshifts
    """
    sigz_c = np.interp(muz_c, z_binc, sigz_z)
    res =  (1./(sigz_c * np.sqrt(2 * np.pi))) * np.exp(-0.5*((z-muz_c)/sigz_c)**2)
    return res

clus_zr(z, r, muz_c, bin_zc, sigz_z, rx)

Evaluate the cluster model at given redshifts and physical radii

Parameters:

Name Type Description Default
z array_like

Array of redshifts

 required
r array_like

Array of physical radii

 required
muz_c float

Mean of the cluster redshift distribution

 required
bin_zc array_like

Array of redshift bins

 required
sigz_z array_like

Array of standard deviation of the galaxy photometric redshift (sigma_zphot(z))

 required
rx float

Scale radius of the cluster model

 required

Returns:

Name Type Description
clus_model array_like

Array of the cluster model evaluated at the input redshifts and physical radii
Source code in photXclus/model.py 
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
@nb.njit()
def clus_zr(z, r, muz_c, bin_zc, sigz_z, rx):
    """
    Evaluate the cluster model at given redshifts and physical radii

    Parameters
    ----------
    z : array_like
        Array of redshifts
    r : array_like
        Array of physical radii
    muz_c : float
        Mean of the cluster redshift distribution
    bin_zc : array_like
        Array of redshift bins
    sigz_z : array_like
        Array of standard deviation of the galaxy photometric redshift (sigma_zphot(z))
    rx : float
        Scale radius of the cluster model

    Returns
    -------
    clus_model : array_like
        Array of the cluster model evaluated at the input redshifts and physical radii
    """
    res = clus_z(z, muz_c, bin_zc, sigz_z) * clus_r(r, rx)
    return res

field_r(r, C_Omega)

Evaluate the field radial model at given physical radii

Parameters:

Name Type Description Default
r array_like

Array of physical radii

 required
C_Omega float

Area of the field of view in degrees^-2

 required

Returns:

Name Type Description
field_model array_like

Array of the field model evaluated at the input physical radii
Source code in photXclus/model.py 
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
@nb.njit()
def field_r(r, C_Omega):
    """
    Evaluate the field radial model at given physical radii

    Parameters
    ----------
    r : array_like
        Array of physical radii
    C_Omega : float
        Area of the field of view in degrees^-2

    Returns
    -------
    field_model : array_like
        Array of the field model evaluated at the input physical radii
    """
    return 2*np.pi*r / C_Omega

field_z(z, binz_ampl_f, bin_z, sigz_z)

Evaluate the field redshift model at given redshifts

Parameters:

Name Type Description Default
z array_like

Array of redshifts

 required
binz_ampl_f array_like

Amplitude of the field model at each redshift bin

 required
bin_z array_like

Array of redshift bins

 required
sigz_z array_like

Array of standard deviation of the galaxy photometric redshift (sigma_zphot(z))

 required

Returns:

Name Type Description
field_model array_like

Array of the field model evaluated at the input redshifts
Source code in photXclus/model.py 
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
@nb.njit()
def field_z(z, binz_ampl_f, bin_z, sigz_z):
    """
    Evaluate the field redshift model at given redshifts

    Parameters
    ----------
    z : array_like
        Array of redshifts
    binz_ampl_f : array_like
        Amplitude of the field model at each redshift bin
    bin_z : array_like
        Array of redshift bins
    sigz_z : array_like
        Array of standard deviation of the galaxy photometric redshift (sigma_zphot(z))

    Returns
    -------
    field_model : array_like
        Array of the field model evaluated at the input redshifts
    """
    iz = np.digitize(z,bin_z)-1
    return binz_ampl_f[iz]

get_pclus_grid(zv, rv, muz_c, bin_zc, sigz_z, rx)

Generates a grid of cluster models evaluated at the input points (z, r).

Parameters:

Name Type Description Default
zv array_like

Array of redshifts.

 required
rv array_like

Array of physical radii.

 required
muz_c float

Mean of the cluster redshift distribution.

 required
bin_zc array_like

bin of the cluster redshift parameter.

 required
sigz_z float

Standard deviation of the cluster redshift distribution.

 required
rx float

X-ray radius of the cluster in units of some X-ray radius estimate.

 required

Returns:

Name Type Description
p array_like

Array of cluster models evaluated at the input (z, r).
Source code in photXclus/likelihood.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
@nb.njit()
def get_pclus_grid(zv, rv, muz_c, bin_zc, sigz_z, rx):
    """
    Generates a grid of cluster models evaluated at the input points (z, r).

    Parameters
    ----------
    zv : array_like
        Array of redshifts.
    rv : array_like
        Array of physical radii.
    muz_c : float
        Mean of the cluster redshift distribution.
    bin_zc : array_like
        bin of the cluster redshift parameter.
    sigz_z : float
        Standard deviation of the cluster redshift distribution.
    rx : float
        X-ray radius of the cluster in units of some X-ray radius estimate.

    Returns
    -------
    p : array_like
        Array of cluster models evaluated at the input (z, r).
    """
    p = np.zeros_like(zv)
    for i in range(len(zv)):
        p[i] = clus_zr(zv[i], rv[i], 
                       muz_c, bin_zc, sigz_z, rx)
    return p

log_likelihood(params, data_in, data_out, other_params, bin_z, sigz_z, zr_grid)

Calculates the log-likelihood of the marked Poisson process associated with a (cluster fov, field fov) = (data_in, data_out) pair of datasets for a given set of parameters params.

Parameters:

Name Type Description Default
params array_like

Array of parameters to be evaluated. The first three elements are the number of galaxies in the cluster (ngal_c), the mean of the cluster redshift distribution (muz_c), and the number of field galaxies (ngal_f). The remaining elements are the amplitudes of the field galaxy bins.

 required
data_in array_like

Array of (z, r) points in the cluster line of sight.

 required
data_out array_like

Array of (z, r) points control line of sight (outside the cluster).

 required
other_params tuple

Tuple containing additional parameters: Omega_in, Omega_out, zmin, zmax, r0, rx, integr_clus.

 required
bin_z array_like

Array of bin edges for the field galaxy bins.

 required
sigz_z float

Standard deviation of the cluster redshift distribution.

 required
zr_grid array_like

Array of (z, r) points in the cluster line of sight for evaluating the integral of the cluster model.

 required

Returns:

Name Type Description
lnL float

The log-likelihood of the given parameters.
Source code in photXclus/likelihood.py 
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
@nb.njit()
def log_likelihood(params, data_in, data_out, other_params, bin_z, sigz_z, zr_grid):
    """
    Calculates the log-likelihood of the marked Poisson process associated 
    with a (cluster fov, field fov) = (data_in, data_out) pair of datasets 
    for a given set of parameters params.

    Parameters
    ----------
    params : array_like
        Array of parameters to be evaluated. The first three elements are the number 
        of galaxies in the cluster (ngal_c), 
        the mean of the cluster redshift distribution (muz_c), 
        and the number of field galaxies (ngal_f). 
        The remaining elements are the amplitudes of the field galaxy bins.
    data_in : array_like
        Array of (z, r) points in the cluster line of sight.
    data_out : array_like
        Array of (z, r) points control line of sight (outside the cluster).
    other_params : tuple
        Tuple containing additional parameters: Omega_in, Omega_out, zmin, zmax, r0, rx, integr_clus.
    bin_z : array_like
        Array of bin edges for the field galaxy bins.
    sigz_z : float
        Standard deviation of the cluster redshift distribution.
    zr_grid : array_like
        Array of (z, r) points in the cluster line of sight for evaluating the integral 
        of the cluster model.

    Returns
    -------
    lnL : float
        The log-likelihood of the given parameters.
    """
    Omega_in, Omega_out, zmin, zmax, r0, rx, integr_clus = other_params

    z_in, r_in = data_in[:, 0], data_in[:, 1]
    z_out, r_out = data_out[:, 0], data_out[:, 1]

    ngal_c, muz_c, ngal_f = params[:3]
    binz_ampl_f = params[3:]

    nbinz=len(bin_z)-1
    bin_zc = np.array([0.5*(bin_z[i]+bin_z[i+1]) for i in range(nbinz)])

    ### 1. normalize PDFs of model.

    ### integral 2D
    zv, rv = zr_grid
    dz = zv[0,1]-zv[0,0]
    dr = rv[1,0]-rv[0,0] 
    L_grid_in = get_pclus_grid(zv, rv, muz_c, bin_zc, sigz_z, rx)
    integ_in =  np.trapz(np.trapz(L_grid_in, dx=dz), dx=dr)

    phi_c = ngal_c / integ_in

    ##1.2 g(z) == field(z)
    integz_field = np.sum(binz_ampl_f) * (bin_zc[1]-bin_zc[0])
    binz_ampl_f /= integz_field

    ### 1.3 full model 
    proba_in = clus_field_zr(z_in, r_in, phi_c, muz_c, bin_zc, sigz_z, 
                                ngal_f, binz_ampl_f, bin_z, Omega_in, rx, 1.)
    proba_out = clus_field_zr(z_out, r_out, phi_c, muz_c, bin_zc, sigz_z, 
                                 ngal_f, binz_ampl_f, bin_z, Omega_out, rx, 0.)
    lnL = np.sum(np.log(proba_in)) + np.sum(np.log(proba_out))

    s_in = (ngal_c + Omega_in * ngal_f) 
    s_out = Omega_out * ngal_f

    lnL -= ( s_in + s_out)

    if np.isnan(lnL):
        print(lnL, s_in, params)

    return lnL

log_prior(theta, params_prior)

Calculates the log-prior probability of the given set of parameters theta.

Parameters:

Name Type Description Default
theta array_like

Array of parameters to be evaluated. The first three elements are the number of galaxies in the cluster (ngal_c), the mean of the cluster redshift distribution (muz_c), and the number of field galaxies (ngal_f). The remaining elements are the amplitudes of the field galaxy bins.

 required
params_prior tuple

Tuple containing additional parameters: ngal_c_max, ngal_f_max, field_bin_size, zrx_min, zrx_max.

 required

Returns:

Name Type Description
lnPrior float

The log-prior probability of the given parameters.
Source code in photXclus/likelihood.py 
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
def log_prior(theta, params_prior):
    """
    Calculates the log-prior probability of the given set of parameters theta.

    Parameters
    ----------
    theta : array_like
        Array of parameters to be evaluated. The first three elements are the number 
        of galaxies in the cluster (ngal_c), the mean of the cluster redshift distribution (muz_c), 
        and the number of field galaxies (ngal_f). 
        The remaining elements are the amplitudes of the field galaxy bins.
    params_prior : tuple
        Tuple containing additional parameters: ngal_c_max, ngal_f_max, 
        field_bin_size, zrx_min, zrx_max.

    Returns
    -------
    lnPrior : float
        The log-prior probability of the given parameters.
    """
    ngal_c, muz_c, ngal_f = theta[:3]
    binz_ampl_f = theta[3:]
    ngal_c_max, ngal_f_max, field_bin_size, zrx_min, zrx_max = params_prior
    if (np.all(0 <= binz_ampl_f) and np.all(binz_ampl_f <= 10000)
        and 1. <= ngal_c <= ngal_c_max 
        and 0. < muz_c < 1.5-field_bin_size
        and 0 <= ngal_f <= ngal_f_max ):
        return logprior_z_rx(muz_c, zrx_min, zrx_max)
    return -np.inf

logprior_z_rx(z, zmin, zmax)

Calculates the log-prior probability of the given redshift z.

Parameters:

Name Type Description Default
z float

Redshift to be evaluated.

 required
zmin float

Minimum allowed redshift.

 required
zmax float

Maximum allowed redshift.

 required

Returns:

Name Type Description
lnPrior float

The log-prior probability of the given redshift.
Source code in photXclus/likelihood.py 
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
def logprior_z_rx(z, zmin, zmax):
    """
    Calculates the log-prior probability of the given redshift z.

    Parameters
    ----------
    z : float
        Redshift to be evaluated.
    zmin : float
        Minimum allowed redshift.
    zmax : float
        Maximum allowed redshift.

    Returns
    -------
    lnPrior : float
        The log-prior probability of the given redshift.
    """
    if (z >= zmin) & (z <= zmax):
        res = 0
    if z > zmax:
        res = -0.5*(z-zmax)**2/((0.025*(1+zmax))**2)
    if z < zmin:
        res = -0.5*(z-zmin)**2/((0.025*(1+zmin))**2)
    return res

plum_r(r, rx)

Evaluate the Plum model at given physical radii

Parameters:

Name Type Description Default
r array_like

Array of physical radii

 required
rx float

Scale radius of the model

 required

Returns:

Name Type Description
plum_model array_like

Array of the Plum model evaluated at the input physical radii
Source code in photXclus/model.py 
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
@nb.njit()
def plum_r(r, rx):
    """
    Evaluate the Plum model at given physical radii

    Parameters
    ----------
    r : array_like
        Array of physical radii
    rx : float
        Scale radius of the model

    Returns
    -------
    plum_model : array_like
        Array of the Plum model evaluated at the input physical radii
    """
    rs = rx/10.
    res =  ( 1 / np.sqrt(1 + (r/rs)**2) ) - (1 / np.sqrt(1 + (10)**2))
    return res

psi_r(r, rs)

Evaluate the dimensionless potential at given physical radii

Parameters:

Name Type Description Default
r array_like

Array of physical radii

 required
rs float

Scale radius of the potential

 required

Returns:

Name Type Description
psi array_like

Array of the dimensionless potential evaluated at the input physical radii
Source code in photXclus/model.py 
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
@nb.njit()
def psi_r(r, rs):    
    """
    Evaluate the dimensionless potential at given physical radii

    Parameters
    ----------
    r : array_like
        Array of physical radii
    rs : float
        Scale radius of the potential

    Returns
    -------
    psi : array_like
        Array of the dimensionless potential evaluated at the input physical radii
    """
    denom = r/rs * (1 + r/rs)**2
    return 1. / denom

psi_r_cut(r, rs)

Evaluate the dimensionless potential at given physical radii, cut off at a certain radius

Parameters:

Name Type Description Default
r array_like

Array of physical radii

 required
rs float

Scale radius of the potential

 required

Returns:

Name Type Description
psi array_like

Array of the dimensionless potential evaluated at the input physical radii, cut off at a certain radius
Source code in photXclus/model.py 
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
@nb.njit()
def psi_r_cut(r, rs):    
    """
    Evaluate the dimensionless potential at given physical radii, cut off at a certain radius

    Parameters
    ----------
    r : array_like
        Array of physical radii
    rs : float
        Scale radius of the potential

    Returns
    -------
    psi : array_like
        Array of the dimensionless potential evaluated at the input physical radii, cut off at a certain radius
    """
    rcut = 0.01*3*rs
    sigma0 = psi_r(rcut, rs) / (2*np.pi*rcut)
    res = psi_r(r, rs)
    i = 0
    while r[i] < rcut:
        res[i] = 2*np.pi*r[i]*sigma0
        i += 1
    return res