dustem_marginalize_grid.pro
6.77 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
PRO dustem_marginalize_grid,fixed_parameters_description $
,fixed_parameters_values $
,seds $
,Nseds $
,params $
,Nparams $
,table_params $
,table_pmins $
,table_pmaxs $
,method=method $
,help=help
;+
; NAME:
; dustem_marginalize_grid
; PURPOSE:
; marginalize a grid on given fixed parameter
; CATEGORY:
; Dustem
; CALLING SEQUENCE:
; dustem_marginalize_grid,fixed_parameters_description,fixed_parameters_values,seds,Nseds,Nfilters,Nparams,table_params,table_pmins,table_pmaxs
; INPUTS:
; fixed_parameters_description : fixed parameters description
; fixed_parameters_values : fixed parameters values
; seds : sed grid to be marginalized (also an output)
; Nseds : Number of seds (also an output)
; params : sed grid parameter values to be marginalized (also an output)
; Nparams : Number of parameters in sed grid (also an output)
; table_params : description of parameters in sed grid (also an output)
; table_pmins : min values of parameters in sed grid (also an output)
; table_pmaxs : max values of parameters in sed grid (also an output)
; OPTIONAL INPUT PARAMETERS:
; NONE
; OUTPUTS:
; None
; OPTIONAL OUTPUT PARAMETERS:
; seds : sed grid to be marginalized (also an output)
; Nseds : Number of seds (also an output)
; Nparams : Number of parameters in sed grid (also an output)
; table_params : description of parameters in sed grid (also an output)
; table_pmins : min values of parameters in sed grid (also an output)
; table_pmaxs : max values of parameters in sed grid (also an output)
; ACCEPTED KEY-WORDS:
; help = If set, print this help
; COMMON BLOCKS:
; None
; SIDE EFFECTS:
; None
; RESTRICTIONS:
; ;Caution: grids should have regularly spaced parameter values for this to work
; PROCEDURE:
; None
; EXAMPLES
;
; MODIFICATION HISTORY:
; Written by J.-Ph. Bernard (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_marginalize_grid'
goto,the_end
ENDIF
N_fixed_params=n_elements(fixed_parameters_description)
;Nparams=!dustem_grid.Nparams
use_method='NEAREST'
IF keyword_set(method) THEN BEGIN
use_method=method
ENDIF
Nfilters=(*!dustem_grid).Nfilters
filters=(tag_names(((*!dustem_grid).st_seds)[0]))[0:Nfilters-1]
FOR i=0L,Nfilters-1 DO BEGIN
filters[i]=strmid(filters[i],1,1000)
ENDFOR
;Nfilters=n_elements(filters)
CASE use_method OF
'NEAREST': BEGIN
FOR i=0L,N_fixed_params-1 DO BEGIN
ind=where(table_params EQ fixed_parameters_description[i],count,complement=indc,ncomplement=nc)
IF count NE 0 THEN BEGIN
;fixed_parameter_value=fixed_parameters_values[ind[0]]
fixed_parameter_value=fixed_parameters_values[i]
;compute distance between table parameter and the requested fixed value
distance=abs(fixed_parameter_value-params[*].(ind[0]))
indd=where(distance EQ min(distance),ccount)
;print,ccount
;2197
;stop
IF ccount EQ 0 THEN BEGIN
message,'This should not happen',/continue
stop
ENDIF
seds=seds[indd]
params=params[indd]
Nseds=ccount
ENDIF
ENDFOR
END
'RELLIKELYHOOD':BEGIN
Nfilters=n_tags(seds[0])-1
FOR i=0L,N_fixed_params-1 DO BEGIN
ind=where(table_params EQ fixed_parameters_description[i],count,complement=indc,ncomplement=nc)
params_v=params[*].(ind[0])
sigma=(params_v*0.1)^2
IF count NE 0 THEN BEGIN
fixed_parameter_value=fixed_parameters_values[i]
print,min(params_v),fixed_parameter_value,max(params_v)
;compute chi2 between table parameter and the requested fixed value
chi2=(fixed_parameter_value-params_v)^2/sigma
min_chi2=min(chi2)
weights=exp(-1.d0*chi2/min_chi2/2.) ;relative Likelyhood
;cgplot,params_v,weights,psym=-4,/ylog
;distance=abs(fixed_parameter_value-params[*].(ind[0]))
;indd=where(distance EQ min(distance),ccount)
indd=where(chi2 EQ min(chi2),ccount)
;print,ccount
;stop
IF ccount EQ 0 THEN BEGIN
message,'This should not happen',/continue
stop
ENDIF
in_seds=seds[indd]
out_seds=in_seds
tot_weights=total(weights)
;stop
;FOR i=0L,count-1 DO BEGIN
ind_not_zero=where(weights NE 0.,count_not_zero)
FOR k=0L,n_elements(out_seds)-1 DO BEGIN
FOR j=0L,Nfilters-1 DO BEGIN
out_seds[k].(j)=total(seds[ind_not_zero].(j)*weights[ind_not_zero])/tot_weights
ENDFOR
ENDFOR
;stop
;==== recompute normalization
FOR k=0L,n_elements(out_seds)-1 DO BEGIN
total_value=0.D0
FOR j=0L,Nfilters-1 DO BEGIN
total_value=total_value+out_seds[k].(j)
ENDFOR
out_seds[k].(Nfilters)=total_value
ENDFOR
seds=out_seds
params=params[indd]
params.(ind[0])=fixed_parameter_value
Nseds=ccount
;stop
ENDIF
ENDFOR
END
'PREDICTOR':BEGIN
;stop
FOR i=0L,N_fixed_params-1 DO BEGIN
ind=where(table_params EQ fixed_parameters_description[i],count,complement=indc,ncomplement=nc)
IF count NE 0 THEN BEGIN
fixed_parameter_value=fixed_parameters_values[i]
;compute distance between table parameter and the requested fixed value
distance=abs(fixed_parameter_value-params[*].(ind[0]))
indd=where(distance EQ min(distance),ccount) ;This will be ne thew number of parameter combinations
;print,ccount
out_seds=seds[indd]
FOR k=0L,ccount-1 DO BEGIN
print,k,ccount
params_v=[params[indd[k]].(0)]
FOR kk=1L,Nparams-1 DO params_v=[params_v,params[indd[k]].(kk)]
params_v[ind[0]]=fixed_parameters_values[i]
parameters_description=(*!dustem_grid).params
parameters_values=params_v
;stop
out_sed=dustem_brute_force_fit_prediction_jpb(parameters_description $
,parameters_values $
;,filters=filters $
,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)
;stop
FOR j=0L,Nfilters-1 DO BEGIN
out_seds[k].(j)=out_sed[j]
ENDFOR
;stop
ENDFOR
seds=out_seds
params=params[indd]
params.(ind[0])=fixed_parameter_value
Nseds=ccount
;stop
ENDIF
ENDFOR
END
ENDCASE
;stop
the_end:
END