FUNCTION dustem_mpfit_run,x,p_min,niter=niter_completed,_EXTRA=_extra ;+ ; NAME: ; dustem_mpfit_run ; ; PURPOSE: ; function used to fit dustem model with dustem_mpfit_sed ; ; CATEGORY: ; Dustem ; ; CALLING SEQUENCE: ; res=dustem_mpfit_run_vg(x,params[,_EXTRA=extra][,cc_sed=cc_sed][,/help]) ; ; INPUTS: ; x = necessary for interface but unused ; params = Dustem Model parameters normalized to initial values ; ; OPTIONAL INPUT PARAMETERS: ; cc_sed = SED color correction factors ; _extra = extra parameters for the plot routine ; ; OUTPUTS: ; None ; ; OPTIONAL OUTPUT PARAMETERS: ; None ; ; ACCEPTED KEY-WORDS: ; help = If set, print this help ; ; COMMON BLOCKS: ; None ; ; SIDE EFFECTS: ; SED and model are plotted if !dustem_show_plot ; ; RESTRICTIONS: ; The dustem idl wrapper must be installed ; ; PROCEDURE: ; None ; ; EXAMPLES ; res=dustem_mpfit_run_vg(x,params[,_EXTRA=extra][,cc_sed=cc_sed][,/help]) ; ; MODIFICATION HISTORY: ; Adapted to extinction and polarization by Vincent Guillet (IAS) ; Written by J.-Ph. Bernard ; Evolution details on the DustEMWrap gitlab. ; See http://dustemwrap.irap.omp.eu/ for FAQ and help. ;- ;note that x is not used. IF keyword_set(help) THEN BEGIN doc_library,'dustem_mpfit_run' sed=0. goto,the_end ENDIF p_dim=p_min*(*(*!dustem_fit).param_init_values) IF !run_ionfrac gt 0. THEN toto = dustem_create_ionfrac() ;RUN THE MODEL AND GET THE SPECTRUM dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st=st ;st is an output... st_model=dustem_read_all(!dustem_soft_dir) chi2_sed = 0 chi2_ext = 0 chi2_qsed =0 chi2_used =0 chi2_qext =0 chi2_uext =0 rchi2_sed = 0 rchi2_ext = 0 rchi2_qsed =0 rchi2_used =0 rchi2_qext =0 rchi2_uext =0 wrchi2_sed = 0 wrchi2_ext = 0 wrchi2_qsed =0 wrich2_used =0 wrchi2_qext =0 wrich2_uext =0 n_sed = 0 n_ext = 0 n_qsed =0 n_used =0 n_qext =0 n_uext =0 win = 0 tagnames=tag_names((*!dustem_data)) dustem_out = [0] ind_data_tag=dustem_data_tag(count=count_data_tag) FOR i_tag=0,count_data_tag-1 DO BEGIN CASE tagnames(ind_data_tag(i_tag)) OF 'SED' : BEGIN dustem_sed = dustem_compute_sed(p_dim,st=st,SED_spec=SED_spec) ;help,*(*!dustem_fit).CURRENT_PARAM_VALUES ;stop chi2_sed = total(((*(*!dustem_data).sed).values-dustem_sed)^2/(*(*!dustem_data).sed).sigma^2) n_sed = n_elements((*(*!dustem_data).sed).values) rchi2_sed = chi2_sed / n_sed wrchi2_sed = rchi2_sed * (*!fit_rchi2_weight).sed message,"chi2 SED = "+strtrim(chi2_sed,2)+" rchi2 SED = "+strtrim(rchi2_sed,2),/continue;," weighted rchi2 SED=",wrchi2_sed ;print,'wl ',(*(*!dustem_data).sed).wav ;print,"Delta SED per wl: ",(*(*!dustem_data).sed).values-dustem_sed ;print,"error per wl: ",(*(*!dustem_data).sed).sigma ;print,"chi2 per wl: ",((*(*!dustem_data).sed).values-dustem_sed)^2/(*(*!dustem_data).sed).sigma^2 alpha = total(((*(*!dustem_data).sed).values*dustem_sed)/(*(*!dustem_data).sed).sigma^2) $ / total(((*(*!dustem_data).sed).values)^2/(*(*!dustem_data).sed).sigma^2) beta = total(((*(*!dustem_data).sed).values*dustem_sed)/(*(*!dustem_data).sed).sigma^2) $ / total(dustem_sed^2/(*(*!dustem_data).sed).sigma^2) print,"Best SED chi2 obtained if f_HI = ",alpha, " or if we multiply model by ",beta ;These two lines were put here so that the chi2 is updated everytime we plot ;Moving it to the bottom so that all chi2 are taken into account ; (*!dustem_fit).chi2=chi2_sed ; (*!dustem_fit).rchi2=rchi2_sed dustem_out = [ dustem_out, dustem_sed ] END 'EXT' : BEGIN dustem_ext = dustem_compute_ext(p_dim,st=st,EXT_spec=EXT_spec) chi2_ext = total(((*(*!dustem_data).ext).values-dustem_ext)^2/(*(*!dustem_data).ext).sigma^2) n_ext = n_elements((*(*!dustem_data).ext).values) rchi2_ext = chi2_ext / n_ext wrchi2_ext = rchi2_ext * (*!fit_rchi2_weight).ext print,"chi2 EXT = ",chi2_ext," rchi2 EXT = ",rchi2_ext;," weighted rchi2 EXT=",wrchi2_ext ;print,"chi2 per wl: ",((*(*!dustem_data).ext).values-dustem_ext)^2/(*(*!dustem_data).ext).sigma^2 dustem_out = [ dustem_out, dustem_ext ] END 'POLEXT': BEGIN ;Haven't investigated !run_circ=1 yet... ;---------------the procedure within this block has not been updated for a vert long time IF !dustem_show_plot NE 0 and !run_circ then begin win+=1 dustem_plot_circ,st,win=win,xr=[0.1,5],aligned=st_model.align.grains.aligned,yr=[-0.03,0.03] endif ;-------------------------------------------------------------------------------------- IF !run_lin then begin dustem_polext = dustem_compute_polext(p_dim,st=st,POLEXT_spec=POLEXT_spec,SPEXT_spec=SPEXT_spec,dustem_fpolext=dustem_fpolext,dustem_ext=dustem_ext) ENDIF END 'QSED': BEGIN stokes = dustem_compute_stokes(p_dim,st=st,Q_spec=Q_spec,U_spec=U_spec,PSI_spec=PSI_spec,dustem_psi_em=dustem_psi_em) dustem_qsed = stokes[0] ;taking care of the color correction for spass ; SHOULD NOT LEAVE THIS ON ;dustem_qsed(-1)=dustem_qsed(-2) ; testspass1 = ((*(*!dustem_data).qsed).filt_names)(-1) eq 'SPASS1' ; if testspass1 then dustem_qsed(-1)=dustem_qsed(-2) ; chi2_qsed = total(((*(*!dustem_data).qsed).values-dustem_qsed)^2/(*(*!dustem_data).qsed).sigma^2) n_qsed = n_elements((*(*!dustem_data).qsed).values) rchi2_qsed = chi2_qsed / n_qsed wrchi2_qsed = rchi2_qsed * (*!fit_rchi2_weight).qsed dustem_out = [ dustem_out, dustem_qsed ] ; classic command message,"chi2 QSED = "+strtrim(chi2_qsed,2)+" rchi2 QSED = "+strtrim(rchi2_qsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed print,"chi2 per wl: ",((*(*!dustem_data).qsed).values-dustem_qsed)^2/(*(*!dustem_data).qsed).sigma^2 END 'USED': BEGIN dustem_used = stokes[1] chi2_used = total(((*(*!dustem_data).used).values-dustem_used)^2/(*(*!dustem_data).used).sigma^2) n_used = n_elements((*(*!dustem_data).used).values) rchi2_used = chi2_used / n_used wrchi2_used = rchi2_used * (*!fit_rchi2_weight).used dustem_out = [ dustem_out, dustem_used] ; classic command message,"chi2 USED = "+strtrim(chi2_used,2)+" rchi2 USED = "+strtrim(rchi2_used,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed print,"chi2 per wl: ",((*(*!dustem_data).used).values-dustem_used)^2/(*(*!dustem_data).used).sigma^2 END 'POLSED': BEGIN IF !run_pol and !run_lin THEN BEGIN ;the keyword for dustem_sed here is used so that the function uses it to compute the polarization fraction instead of ;running dustem_compute_sed again in the procedure. Which it used to do. dustem_polsed = dustem_compute_polsed(p_dim,st=st,P_spec=P_spec,SP_spec=SP_spec,dustem_polfrac=dustem_polfrac,dustem_sed=dustem_sed) ENDIF END 'QEXT': BEGIN ;UPDATED Stokext = dustem_compute_stokext(p_dim,st=st,QEXT_spec=QEXT_spec,UEXT_spec=UEXT_spec,PSIEXT_spec=PSIEXT_spec,dustem_psi_ext=dustem_psi_ext) dustem_qext = stokext[0] chi2_qext =total(((*(*!dustem_data).qext).values-dustem_qext)^2/(*(*!dustem_data).qext).sigma^2) n_qext = n_elements((*(*!dustem_data).qext).values) rchi2_qext = chi2_qext / n_qext wrchi2_qext = rchi2_qext * (*!fit_rchi2_weight).qext dustem_out = [ dustem_out, dustem_qext] ; classic command message,"chi2 QEXT = "+strtrim(chi2_qext,2)+" rchi2 QEXT = "+strtrim(rchi2_qext,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed print,"chi2 per wl: ",((*(*!dustem_data).qext).values-dustem_qext)^2/(*(*!dustem_data).qext).sigma^2 END ; 'UEXT': BEGIN ;UPDATED dustem_uext = stokext[1] chi2_uext =total(((*(*!dustem_data).uext).values-dustem_uext)^2/(*(*!dustem_data).uext).sigma^2) n_uext = n_elements((*(*!dustem_data).uext).values) rchi2_uext = chi2_uext / n_uext wrchi2_uext = rchi2_uext * (*!fit_rchi2_weight).uext dustem_out = [ dustem_out, dustem_uext] ; classic command message,"chi2 UEXT = "+strtrim(chi2_uext,2)+" rchi2 UEXT = "+strtrim(rchi2_uext,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed print,"chi2 per wl: ",((*(*!dustem_data).uext).values-dustem_uext)^2/(*(*!dustem_data).uext).sigma^2 END ELSE : ;Do nothing ENDCASE ENDFOR ;=== computing chi^2 DOF = n_elements(p_min) total_chi2 = 0 total_rchi2 = 0 total_n = 0 if (*!fit_rchi2_weight).ext ne 0 then begin total_n += n_ext total_chi2 += chi2_ext total_rchi2 += rchi2_ext endif if (*!fit_rchi2_weight).sed ne 0 then begin total_n += n_sed total_chi2 += chi2_sed total_rchi2 += rchi2_sed endif if !run_pol then begin if tag_exist(*!fit_rchi2_weight,'qsed') then begin ;takes into account the activation$ ;of the !run_pol keyword without the user knowing it. ;ie: the use of a model with polarization and without polarization data or the user wanting to fit polarization. if (*!fit_rchi2_weight).qsed ne 0 then begin total_n += n_qsed total_chi2 += chi2_qsed total_rchi2 += rchi2_qsed endif if (*!fit_rchi2_weight).used ne 0 then begin total_n += n_used total_chi2 += chi2_used total_rchi2 += rchi2_used endif endif if tag_exist(*!fit_rchi2_weight,'qext') then begin if (*!fit_rchi2_weight).qext ne 0 then begin total_n += n_qext total_chi2 += chi2_qext total_rchi2 += rchi2_qext endif if (*!fit_rchi2_weight).uext ne 0 then begin total_n += n_uext total_chi2 += chi2_uext total_rchi2 += rchi2_uext endif endif endif total_wchi2 = chi2_ext * (*!fit_rchi2_weight).ext $ + chi2_sed * (*!fit_rchi2_weight).sed if !run_pol then begin if tag_exist(*!fit_rchi2_weight,'qsed') THEN total_wchi2+=chi2_qsed * (*!fit_rchi2_weight).qsed $ +chi2_used * (*!fit_rchi2_weight).used if tag_exist(*!fit_rchi2_weight,'qext') THEN total_wchi2+=chi2_qext * (*!fit_rchi2_weight).qext $ +chi2_used * (*!fit_rchi2_weight).used ; total_wchi2 += $;chi2_polext * !fit_rchi2_weight.polext $ ; ; + chi2_polsed * !fit_rchi2_weight.polsed $ ; ; + chi2_polfrac * !fit_rchi2_weight.polfrac $ ; + chi2_qsed * (*!fit_rchi2_weight).qsed $ ; + chi2_used * (*!fit_rchi2_weight).used $ ; + chi2_qext * (*!fit_rchi2_weight).qext $ ; + chi2_uext * (*!fit_rchi2_weight).uext endif ;Using this expression for now for the weighting. (*!dustem_fit).chi2 = total_wchi2 (*!dustem_fit).rchi2 = total_wchi2/(total_n-DOF) message,'======================================================================',/continue print,"Total weighted chi2 driving mpfitfun = ", total_wchi2 print,"Nb of data points = ", total_n print,"DOF = ", DOF print,"Reduced chi2 EXT/SED/POLEXT/POLSED/POLFRAC",rchi2_ext,rchi2_sed;,rchi2_polext;,rchi2_polsed,rchi2_polfrac print,"Weighted reduced chi2 EXT/SED/POLEXT/POLSED/POLFRAC",wrchi2_ext,wrchi2_sed;,wrchi2_polext;,wrchi2_polsed,wrchi2_polfrac if !run_pol then print,"Weights",(*!fit_rchi2_weight).ext,(*!fit_rchi2_weight).sed $;,!fit_rchi2_weight.polext $;,!fit_rchi2_weight.polsed,!fit_rchi2_weight.polfrac $ else print,"Weights",(*!fit_rchi2_weight).ext,(*!fit_rchi2_weight).sed print,"Mean weighted reduced chi2 = ", total_wchi2/(total_n-DOF) print,"Unweighted reduced Total chi2 = ", total_chi2/(total_n-DOF) print,"Unweighted mean reduced chi2 = ", total_rchi2/4*total_n/(total_n-DOF);total_rchi2/4*total_n/(total_n-DOF) original ; I do not understand this line. Maybe redo some computation. message,'======================================================================',/continue ;the following is needed under fawlty (not IDL) which message,/reset does not work like in IDL !Error_State.msg='' If !dustem_show_plot THEN BEGIN ;=== plotting the fit if !dustem_iter.act NE !dustem_iter.prv or !dustem_iter.prv EQ 1 then begin IF !dustem_noobj EQ 0 THEN BEGIN dustemwrap_plot, p_dim, $ st = st, $ dustem_sed = dustem_sed, $ SED_spec = SED_spec, $ dustem_qsed = dustem_qsed, $ Q_spec = Q_spec, $ dustem_used = dustem_used, $ U_spec = U_spec, $ dustem_polsed = dustem_polsed, $ P_spec = P_spec, $ dustem_polfrac = dustem_polfrac, $ SP_spec = SP_spec, $ dustem_psi_em = dustem_psi_em, $ PSI_spec = PSI_spec, $ dustem_ext = dustem_ext, $ EXT_spec = EXT_spec, $ dustem_qext = dustem_qext, $ QEXT_spec = QEXT_spec, $ dustem_uext = dustem_uext, $ UEXT_spec = UEXT_spec, $ dustem_polext = dustem_polext, $ POLEXT_spec = POLEXT_spec, $ dustem_fpolext = dustem_fpolext, $ SPEXT_spec = SPEXT_spec, $ dustem_psi_ext = dustem_psi_ext, $ PSIEXT_spec = PSIEXT_spec, $ _extra=_extra ENDIF ELSE BEGIN dustemwrap_plot_noobj, p_dim,$ st = st,$ dustem_sed = dustem_sed, $ SED_spec = SED_spec, $ dustem_qsed = dustem_qsed, $ Q_spec = Q_spec, $ dustem_used = dustem_used, $ U_spec = U_spec, $ dustem_polsed = dustem_polsed, $ P_spec = P_spec, $ dustem_polfrac = dustem_polfrac, $ SP_spec = SP_spec, $ dustem_psi_em = dustem_psi_em, $ PSI_spec = PSI_spec, $ dustem_ext = dustem_ext, $ EXT_spec = EXT_spec, $ dustem_qext = dustem_qext, $ QEXT_spec = QEXT_spec, $ dustem_uext = dustem_uext, $ UEXT_spec = UEXT_spec, $ dustem_polext = dustem_polext, $ POLEXT_spec = POLEXT_spec, $ dustem_fpolext = dustem_fpolext, $ SPEXT_spec = SPEXT_spec, $ dustem_psi_ext = dustem_psi_ext, $ PSIEXT_spec = PSIEXT_spec, $ _extra = _extra ENDELSE if !dustem_iter.act EQ 1 then begin !dustem_iter.prv=0 !dustem_iter.act=0 endif else !dustem_iter.prv=!dustem_iter.act endif ENDIF ;print,"Total reduced chi2 = ", (wrchi2_sed+wrchi2_ext+wrchi2_polext+wrchi2_polsed+wrchi2_polfrac)/(n_sed+n_ext+n_polext+n_polsed+n_polfrac-DOF) ; Reduced chi2 is like if there were the same nb of data points for each curve ;print,"Mean reduced chi2 = ", (wrchi2_sed+wrchi2_ext+wrchi2_polext+wrchi2_polsed+wrchi2_polfrac) * total_n/(total_n-DOF) ; Plot size distribution win+=1 xr=[0.3,5e3] yr=[0.1,100] ;IF !dustem_show_plot NE 0 THEN dustem_plot_sdist,xr=xr,yr=yr,ps=filename,win=win;,/donotclose ; Suppress the artificial [0] at the beginning dustem_out = dustem_out[1:*] ;IF defined(extra) THEN BEGIN ; IF keyword_set(extra.plot_it) THEN BEGIN ;; message,'params='+strtrim(p_min,2),/info,/continue ;; stop ; print,'dustem_mpfit_run_vg: Current parameters values:',p_dim ; ENDIF ;ENDIF the_end: RETURN,dustem_out END