-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfitting.py
262 lines (227 loc) · 7.66 KB
/
fitting.py
1
2
3
4
5
6
7
8
9
10
11
12
13
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
44
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
121
122
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
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
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
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
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy import integrate
import determine_radius as dr
import splashback as sp
plt.style.use("mnras.mplstyle")
H0 = 67.77 * 1000 / 3.09e22
G = 6.67e-11
rho_crit = 3 * H0**2 / (8 * np.pi * G) #kg/m^3
unit_converter = 1.99e30 / (3.09e22**3)
rho_crit = rho_crit / unit_converter
def log_grad_model(log_radii, rho_s, r_s, r_t, alpha, beta, gamma, b_e, S_e):
"""
Density model from DK14.
Parameters
----------
log_radii : array
Log radii of profiles, given in r/R200m.
rho_s : str
Fitting parameter.
r_s : str
Fitting parameter.
r_t : str
Fitting parameter.
alpha : str
Fitting parameter.
beta : str
Fitting parameter.
gamma : str
Fitting parameter.
b_e : str
Fitting parameter.
S_e : str
Fitting parameter.
Returns
-------
log_grad : array
Log-gradient of profiles.
"""
rho_s = 10**rho_s
radii = 10**log_radii
rho_m = 0.307 * rho_crit
rho_inner = rho_s * np.exp(-2/alpha * ((radii/r_s)**alpha - 1))
f_trans = (1 + (radii/r_t)**beta)**(-gamma/beta)
rho_outer = rho_m * (b_e * (radii/5)**(-1*S_e) + 1)
rho_total = rho_inner * f_trans + rho_outer
log_grad = np.gradient(np.log10(rho_total), log_radii)
return log_grad
def density_model(log_radii, rho_s, r_s, r_t, alpha, beta, gamma, b_e, S_e):
"""
DK14 density model.
log_radii : array
Log radii of profiles, given in r/R200m.
rho_s : str
Fitting parameter.
r_s : str
Fitting parameter.
r_t : str
Fitting parameter.
alpha : str
Fitting parameter.
beta : str
Fitting parameter.
gamma : str
Fitting parameter.
b_e : str
Fitting parameter.
S_e : str
Fitting parameter.
Returns
-------
log_density : array
Density from parameters.
"""
rho_s = 10**rho_s
radii = 10**log_radii
rho_m = 0.307 * rho_crit
rho_inner = rho_s * np.exp(-2/alpha * ((radii/r_s)**alpha - 1))
f_trans = (1 + (radii/r_t)**beta)**(-gamma/beta)
rho_outer = rho_m * (b_e * (radii/5)**(-1*S_e) + 1)
rho_total = rho_inner * f_trans + rho_outer
log_density = np.log10(rho_total)
return log_density
def fit_log_models(flm, split, bootstrap=False):
"""
Fits models to
Parameters
----------
flm : TYPE
DESCRIPTION.
split : TYPE
DESCRIPTION.
bootstrap : TYPE, optional
DESCRIPTION. The default is False.
Returns
-------
params : TYPE
DESCRIPTION.
"""
profile = getattr(flm, split+"_profile_DM")
log_profile = getattr(flm, split+"_log_DM")
if bootstrap:
log_profile = getattr(flm, split+"_log_DM_temp")
profile = getattr(flm, split+"_profile_DM_temp")
N_bins = np.shape(log_profile)[0]
guess = np.array([1, 0.15, 1.5, 0.06, 6, 84, 1.5, 0.9], dtype=np.float32)
params = np.zeros((N_bins, 8))
for i in range(N_bins):
guess[0] = np.log10(profile[i,0]*100)
try:
popt, pcov = curve_fit(log_grad_model, np.log10(flm.rad_mid),
log_profile[i,:],
guess, maxfev=10000, method='trf')
params[i,:] = popt
except RuntimeError:
params[i,:] = np.nan
return params
def project_model(radii, params):
"""
Uses density profile parameters to get a projected density profile.
Parameters
----------
radii : array, float
Radii to find profile for.
params : array, float
Fitted parameters for profile
Returns
-------
projected_density : array, float
Projected density profile
"""
R_max = 5
N_rad = len(radii)
N_bins = np.shape(params)[0]
projected_density = np.zeros((N_bins, N_rad))
for i in range(N_bins):
popt = params[i,:]
for j in range(N_rad):
R = radii[j]
integrand = lambda r: 10**(density_model(np.log10(r), popt[0],
popt[1], popt[2], popt[3],
popt[4], popt[5], popt[6],
popt[7])) * r / np.sqrt(r**2 - R**2)
projected_density[i,j], err = integrate.quad(integrand, R, R_max,
epsrel=1e-40)
return projected_density
def find_sort_R(flm, radii, array, names, plot=False):
"""
Finds radii of minima locations in profile. Automatically sorts between
second caustics and splashback features
Parameters
----------
radii : numpy array, float
Array of radii corresponding to profile.
array : numpy array, float
Profile values corresponding to radii locations.
names : list
List of length two. First entry giving the name of the profile, e.g.
DM, gas or SZ. Second entry gives the name of the stacking criteria.
Returns
-------
None.
"""
R_sp, second = dr.depth_cut(radii, array,
cut=-1,
second_caustic="y",
plot=plot)
second_mask = np.where(np.isfinite(second))[0]
for i in range(len(second_mask)):
index = second_mask[i]
if R_sp[index] < second[index]:
larger = second[index]
#smaller = R_sp[index]
R_sp[index] = larger
#second[index] = smaller
return R_sp
def bootstrap_errors(data, split):
"""
Calculates sampling error of splashback radius from fitted and projected
dark matter density profiles.
Parameters
----------
data : obj
simulation data
split : str
Name of criteria used to stack the profiles.
Returns
-------
Rsp_error : array, float
Error values for splashback radius from projected density profiles.
"""
stacking_data = data.DM_density_3D
if split == "mass":
split_data = np.log10(data.M200m)
else:
split_data = getattr(data, split)
split_bins = getattr(data, split+"_bins")
N_bootstrap = 100
not_nan = np.where(np.isfinite(split_data)==True)[0]
bins_sort = np.digitize(split_data[not_nan], split_bins)
N_bins = len(split_bins) - 1
Rsp_error = np.zeros(N_bins)
# depth_error = np.zeros(N_bins)
for i in range(N_bins):
bin_mask = np.where(bins_sort == i+1)[0]
stacked_data = np.zeros((N_bootstrap, data.N_rad))
log_sample = np.zeros((N_bootstrap, data.N_rad))
print(i, len(bin_mask))
for j in range(N_bootstrap):
# Select random sample from bin with replacement
sample = np.random.choice(bin_mask,
size=len(bin_mask),
replace=True)
stacked_data[j,:] = sp.stack_data(stacking_data[not_nan[sample],:])
log_sample = sp.log_gradients(data.rad_mid, stacked_data)
setattr(data, split+"_profile_DM_temp", stacked_data)
setattr(data, split+"_log_DM_temp", log_sample)
params = fit_log_models(data, split, bootstrap=True)
projected_density = project_model(data.rad_mid, params)
projected_model_log_DM = sp.log_gradients(data.rad_mid, projected_density,
smooth=False)
# #Projected splashback model from 3D
R_model = find_sort_R(data, data.rad_mid, projected_model_log_DM,
["model", split])
Rsp_error[i] = np.nanstd(R_model)
return Rsp_error