dustem_brute_force_fit_prediction.pro
8.14 KB
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
FUNCTION dustem_brute_force_fit_prediction,parameters_description,parameters_values,table_name $
,filters $
,reset=reset $
,weighted_seds=weighted_seds $
,best_sed =best_sed $
,best_parameters =best_parameters $
,chi2_p=chi2_p $
,rchi2_p=rchi2_p $
,total_chi2s_p=total_chi2s_p $
,status=status $
,show_sed=show_sed $
,nostop=nostop
;+
; NAME:
; dustem_brute_force_fit_prediction
; PURPOSE:
; Do the prediction of sed for a given input parameters
; CATEGORY:
; DustEM
; CALLING SEQUENCE:
; res=dustem_brute_force_fit(sed,table_name,filters[,/normalize][,fact=][,chi2=][,rchi2=][rchi2][best_sed=][,/show_sed][,/nostop][,/reset])
; INPUTS:
; table_name = name of the fits file to be used as a grid (see dustem_make_sed_table.pro for creation)
; filters = filters to be used in the grid table (only used if !dustem_grid does not exist or if /reset)
; fixed_parameters_description = fixed parameters of the table that should be considered in order to make the prediction of sed
; fixed_parameters_values = fixed parameters values of the table that should be considered in order to make the prediction of sed
; OPTIONAL INPUT PARAMETERS:
; fixed_parameters_description = if set, describes parameters of the table that should be considered fixed.
; fixed_parameters_values = values of parameters of the table that should be considered fixed.
; OUTPUTS:
; res = best values of the free parameters. Not that linear parameters to the sed should be multiplied by fact.
; OPTIONAL OUTPUT PARAMETERS:
; chi2_p = chi2 of the best fit of input parameters
; rchi2_p = rchi2 is the reduce chi2 of the best fit of input parameters
; best_sed = best fit SED matching input parameters
; weighted_seds = likelyhood (1/chi2) weighted best seds values
; total_chi2s_p = all chi2s of fit of input parameters
; ACCEPTED KEY-WORDS:
; help = if set, print this help
; COMMON BLOCKS:
; None
; SIDE EFFECTS:
; None
; RESTRICTIONS:
; The DustEMWrap IDL code must be installed
; PROCEDURE:
; Scans the grid SEDs to find the best fit (lowest chi2)
; EXAMPLES
;
; MODIFICATION HISTORY:
; Written by RM May 2024
; Evolution details on the DustEMWrap gitlab.
; See http://dustemwrap.irap.omp.eu/ for FAQ and help.
;-
IF keyword_set(help) THEN BEGIN
doc_library,'dustem_brute_force_fit_prediction'
params=0.
goto,the_end
ENDIF
;==== test for the existence of a loaded grid and initialize if needed.
;==== only filters contained in the filters variable are kept in the grid, and the sum of grid seds are recomputed accordingly
defsysv,'!dustem_grid',0,exist=exist
IF exist EQ 0 or keyword_set(reset)THEN BEGIN
dustem_define_grid,table_name,filters=filters,status=status
IF status NE 1 THEN BEGIN
params=-1
goto,the_end
ENDIF
ENDIF
;seds=!dustem_grid.seds
seds=(*!dustem_grid).st_seds
params=(*!dustem_grid).st_params
totals=(*!dustem_grid).sed_totals
Nseds=(*!dustem_grid).Nsed
Nfilters=(*!dustem_grid).Nfilters
Nparams=(*!dustem_grid).Nparams
table_params=(*!dustem_grid).params
table_pmins=(*!dustem_grid).pmin_values
table_pmaxs=(*!dustem_grid).pmax_values
wave = dustem_filter2wav(filters)
;=== These are the full arrays needed for calculations
chi2s_p=dblarr(Nseds,Nparams)
grid_seds=dblarr(Nseds,Nfilters)
N_fixed_params=n_elements(parameters_description)
grid_param_values=dblarr(Nseds,Nparams)
err=double(parameters_values*1e-1)
;err = 0.01/100. ; Arbitrary choice of sigma in order to compute the chi2
this_var = la_power(err,2)
; this_var[0]=this_var[0]*1.d-3
;=== compute sed grid Array
FOR i=0L,Nseds-1 DO BEGIN
FOR j=0L,Nfilters-1 DO BEGIN
;grid_seds[i,j]=seds[i].(j+Nparams)
grid_seds[i,j]=seds[i].(j)
ENDFOR
ENDFOR
;=== compute the parameter values Array
FOR i=0L,Nseds-1 DO BEGIN
FOR j=0L,Nparams-1 DO BEGIN
;grid_param_values[i,j]=seds[i].(j)
grid_param_values[i,j]=params[i].(j)
ENDFOR
ENDFOR
marker= 0 ;create a marker in order call just one time dustem_init in order to have default parameters
FOR i=0L,Nparams-1 DO BEGIN
ind=where(parameters_description EQ table_params[i],count,complement=indc,ncomplement=nc)
IF count NE 0 THEN BEGIN
p=params[*].(ind[0])
pval = parameters_values[ind[0]]
chi2s_p[*,ind[0]] = (pval-p)^2/this_var[ind[0]]
ENDIF
IF count EQ 0 THEN BEGIN
marker=marker+1
IF marker EQ 1 THEN BEGIN
dustem_init, model='DBP90'
ENDIF
p=params[*].(i)
pval = dustem_get_param_values_default(/verbose,[table_params[i]])
pval = pval[0]
chi2s_p[*,i] = (pval-p)^2/this_var[0]
ENDIF
ENDFOR
grid_params_bidon = grid_param_values
grid_seds_bidon = grid_seds
all_chi2s_p=chi2s_p
total_chi2s_p = dblarr(Nseds)
total_chi2s_G0= dblarr(Nseds)
FOR i=0L,Nparams-1 DO BEGIN
total_chi2s_p = total_chi2s_p + all_chi2s_p[*,i]
ENDFOR
chi2_p=min(total_chi2s_p)
ind=where(total_chi2s_p EQ chi2_p,count)
grid_params_bidon= grid_params_bidon[ind[0],*]
grid_seds_bidon = grid_seds_bidon[ind[0],*]
chi2s_p = chi2s_p[ind[0],*]
Nfreedom=Nfilters-Nparams-1 ;This is the degree of freedom
rchi2_p=chi2_p/Nfreedom ;This is the reduced rchi2_p
best_parameters = reform(grid_params_bidon) ;This is the best parameters fround in the grid.
best_sed=reform(grid_seds_bidon) ;This is the best SED fround in the grid.
;=== compute Likelyhood weighted seds values
weighted_seds=dblarr(Nfilters)
weights=exp(-total_chi2s_p/2)
tot_weights=total(weights)
IF tot_weights EQ 0. THEN BEGIN
FOR i=0L,5-1 DO BEGIN
weighted_seds[i]=best_sed[i]
ENDFOR
message,'weigthed divergence of sed of PAHs and VSGs! ',/continue
ENDIF ELSE BEGIN
FOR i=0L,5-1 DO BEGIN
weighted_seds[i]=total(reform(1.d0*grid_seds[*,i])*weights)/tot_weights
ENDFOR
ENDELSE
FOR i=1L,Nparams-1 DO BEGIN
total_chi2s_G0 = total_chi2s_G0 + all_chi2s_p[*,i]
ENDFOR
chi2_G0=min( total_chi2s_G0)
ind1=where( total_chi2s_G0 EQ chi2_G0,count)
weights=exp(-all_chi2s_p[ind1,0]/2)
tot_weights=total(weights)
IF tot_weights EQ 0. THEN BEGIN
FOR i=5L,Nfilters-1 DO BEGIN
weighted_seds[i]=best_sed[i]
ENDFOR
message,'weigthed divergence of sed of BGs ! ',/continue
ENDIF ELSE BEGIN
FOR i=5L,Nfilters-1 DO BEGIN
weighted_seds[i]=total(reform(1.d0*grid_seds[ind1,i])*weights)/tot_weights
ENDFOR
ENDELSE
IF keyword_set(show_sed) THEN BEGIN
WINDOW,0,XSIZE=800,YSIZE=800
xtit='Wavelength [mic]'
ytit='SED (weighted_sed=red,best_sed=blue)'
;should also plot uncertainties
minv=la_min([weighted_seds,best_sed])
maxv=la_max([weighted_seds,best_sed])
cgplot,wave,weighted_seds,psym='Filled Circle',color='blue',/ylog,/ysty,xtit=xtit,ytit=ytit,/xlog,yrange=[minv,maxv]
cgoplot,wave+0.5,best_sed,psym='Filled Square',color='red'
; WINDOW,1,XSIZE=800,YSIZE=800
; cgplot, grid_param_values[*,0], 1/all_chi2s_p[*,0],psym='Filled Circle',/xlog,/ylog,ytit='exp(-chi2s/2)',xtit='G0s'
; WINDOW,2,XSIZE=800,YSIZE=800
; cgplot,grid_param_values[*,0], grid_seds[*,3],psym='Filled Circle',/xlog,/ylog,ytit='F1130',xtit='G0s'
; WINDOW,3,XSIZE=800,YSIZE=800
; cgplot,grid_param_values[*,0], grid_seds[*,5],psym='Filled Circle',/xlog,/ylog,ytit='PACS70',xtit='G0s'
; WINDOW,4,XSIZE=800,YSIZE=800
; cgplot,grid_param_values[*,0], DERIV(grid_param_values[*,0],grid_seds[*,3]),psym='Filled Circle',/xlog,/ylog,ytit='D(F1130)/DG0',xtit='G0s'
; WINDOW,5,XSIZE=800,YSIZE=800
; cgplot,grid_param_values[*,0], DERIV(grid_param_values[*,0],grid_seds[*,5]),psym='Filled Circle',/xlog,/ylog,ytit='D(PACS70)/DG0',xtit='G0s'
; WINDOW,6,XSIZE=800,YSIZE=800
; cgplot,grid_seds[*,5], 1/total_chi2s_p,psym='Filled Circle',/xlog,/ylog,ytit='1/chi2s',xtit='PACS70'
; WINDOW,7,XSIZE=800,YSIZE=800
; cgplot,grid_seds[*,3], 1/all_chi2s_p[ind1,0],psym='Filled Circle',/xlog,/ylog,ytit='1/chi2s',xtit='F1130'
; WINDOW,8,XSIZE=800,YSIZE=800
; cgplot,grid_seds[*,2], 1/all_chi2s_p[ind1,0],psym='Filled Circle',/xlog,/ylog,ytit='1/chi2s',xtit='F1000'
;
if not keyword_set(nostop) THEN stop
ENDIF
the_end:
RETURN,weighted_seds
END