dustem_mpfit_data.pro
4.93 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
FUNCTION dustem_mpfit_data,tol=tol,Nitermax=Nitermax,gtol=gtol,function_name=function_name,xtol=xtol,_extra=_extra
;+
; NAME:
; dustem_mpfit_data
; PURPOSE:
; call MPFIT to fit observational datas using ponderated reduced chi2
; between extinction, emission, and polarisation as defined in !fit_rchi2_weight
; CATEGORY:
; Dustem
; CALLING SEQUENCE:
; res=dustem_mpfit_data(tol=tol,Nitermax=Nitermax,gtol=gtol)
; INPUTS:
; None
; OPTIONAL INPUT PARAMETERS:
; OUTPUTS:
; None
; OPTIONAL OUTPUT PARAMETERS:
; ACCEPTED KEY-WORDS:
; help = If set, print this help
; COMMON BLOCKS:
; None
; SIDE EFFECTS:
; None
; RESTRICTIONS:
; The dustem idl wrapper must be installed
; PROCEDURE:
; None
; EXAMPLES
;
; MODIFICATION HISTORY:
; Written by V. Guillet (2012)
; see evolution details on the dustem cvs maintained at CESR
; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-
IF keyword_set(mode) THEN BEGIN
use_mode=strupcase(mode)
ENDIF ELSE BEGIN
use_mode='COMPIEGNE_ETAL2010' ;Default is last dustem model
ENDELSE
disp_str='====================================================='
IF not keyword_set(function_name) THEN BEGIN
func_name='dustem_mpfit_run'
ENDIF ELSE BEGIN
func_name=function_name
ENDELSE
;functargs={plot_it:1}
use_tol=1.e-14 ;Default tolerance value
use_Nitermax=3 ;Default maximum number of iterations
IF keyword_set(tol) THEN use_tol=tol
IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax
; Build array concatenating observational data
if !run_pol then begin
ind_data_tag=dustem_data_tag(count=count_data_tag)
tagnames=tag_names(*!dustem_data)
testpol = tagnames EQ 'POLSED' or tagnames EQ 'POLEXT' or tagnames EQ 'POLFRAC' or tagnames EQ 'PSI_EM' or tagnames EQ 'PSI_EXT' or tagnames EQ 'FPOLEXT'
data_ind=where(~testpol,ctestpol)
if ctestpol ne 0 then begin
match,ind_data_tag,data_ind,ind1,ind2 ;not using match2 because when matching indices, if the first is 0 it doesn't count which is wrong.
;I also only need the points that match.
ind_data_tag = ind_data_tag[ind1]
count_data_tag = n_elements(ind_data_tag)
endif
endif else ind_data_tag=dustem_data_tag(count=count_data_tag)
fields=['wav','values','sigma']
FOR j=0,n_elements(fields)-1 DO BEGIN
instruc = fields(j)+' = ([0' ; we add 0 first to allow for ',' after it
for i=0, count_data_tag-1 do begin $ ; Add each data set
; Adapt sigma to the reduced chi2 ponderation defined in !fit_rchi2_weight between sed, ext, polext and polsed
; Total_reduced_chi2_driving_mpfit = sum_i ( chi2_i * (weight_i/n_i) )
; n : the total number of observational variables in the data
; n_i : the number of obs. variable in data i
; w_i : the weight for data i from !fit_rchi2_weight for the ponderation of reduced chi2
; ex : w_sed = 0.8, w_ext = 0.2 , i.e. rchi2_sed counts 4 times as much as rchi2_ext
; sed : 314 pts, ext : 52 pts, i.e. 6 times much more constraints in sed than in ext
; sigma_sed *= sqrt(314./366/0.8) = 1.03
; sigma_ext *= sqrt(52./366./0.2) = 0.84
; By increasing sigma by a factor w, we diminish the weight of the corresponding point by a factor w^2
if count_data_tag gt 1 and j eq 2 then begin ; sigma ponderation = sqrt(ni/wi)
if !fit_rchi2_weight.(i) gt 0 then begin
weight = sqrt(1./!fit_rchi2_weight.(i) )
endif else begin
weight = 1.d99
endelse
endif else begin
weight = 1.
endelse
;print,"WEIGHT",i,j,weight,n_elements((*!dustem_data.(i)).sigma),"/",n_elements(wav)
instruc += ', ((*(*!dustem_data).(' + strtrim(ind_data_tag(i),2) + ')).'+fields(j)+')*'+strtrim(weight,2)
;print,"INSTRUC",instruc
endfor
instruc += '])(1:*)' ; from 1 to exclude first 0
;print,'-> INSTRUC: ',instruc
toto=execute(instruc) ; Variables wav, values & sigma are created
ENDFOR
nocatch=!dustem_nocatch
p_min = dustem_mpfitfun(func_name,wav,values,sigma, $
parinfo=(*!dustem_parinfo), perror=perror, yfit=yfit, $
bestnorm=bestnorm, dof=dof,quiet=quiet, status=status, xtol=xtol, niter=niter_completed, ftol=use_tol, $
maxiter=use_Nitermax,gtol=gtol,errmsg=errmsg,functargs=_extra,nocatch=nocatch)
message,errmsg,/info
message,'mpfitfun Status='+strtrim(status,2),/info
p_dim=p_min*(*(*!dustem_fit).param_init_values)
p_er_dim=perror*(*(*!dustem_fit).param_init_values)
;PRINT RESULTS
;print,'<°)))>< <°)))>< <°)))>< <°)))>< <°)))>< <°)))>< <°)))>< <°)))><'
print, 'Final parameters : ', p_dim
;=== store best fit values in syst variables
(*!dustem_fit).chi2=bestnorm
(*!dustem_fit).rchi2=bestnorm/dof
(*!dustem_fit).current_param_values=ptr_new(p_dim)
;The above was removed because it seems it was not used anywhere.
;!dustem_fit_sed=ptr_new(dustem_compute_sed(p_dim,st,cont=cont,_extra=extra))
(*!dustem_fit).current_param_errors=ptr_new(p_er_dim)
dustemwrap_plot,p_dim,stp,xr=_extra.xr,/xstyle,yr=_extra.yr,/ysty,/ylog,/xlog,title=_extra.title
RETURN, p_dim
END