dustem_mpfit_data.pro
4.23 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
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 datas
; Order : that of !dustem_data (sed, ext, polext, polsed)
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
;stop
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)
RETURN, p_dim
END