;+ ; PROJECT: RHESSI ; ; NAME: spex_chi2_map_results ; ; PURPOSE: OSPEX routine to analyze the results from the spex::chi2_map method to create a save file ; containing a mapping of chisqr for a sequence of parameter values. Calculates the 1- and 2-sigma ; uncertainties. If plot is requested, the chi-square at each fit is plotted vs the selected fit parameter ; (usually looks like a smooth inverted gaussian), and the parameter values at the chisqr min plus 1 and ; the chisqr min plus 4 give the 1- and 2-sigma error estimates. ; ; EXPLANATION: This is called by the spex::chi2_map method, or can be called directly by the user after ; having created a save file in spex::chi2_map. Note: if the chi-square value doesn't reach the min plus 1 ; or plus 4 level, the min or max of the parameter values is used, and the invalid flag is set (s1_valid, ; s2_valid in out_struct. ; ; CALLING SEQUENCE: ; spex_chi2_map_results ; ; INPUT KEYWORDS: ; savefile - name(s) of savefile(s) to read (written by spex::chi2_map) - scalar or array ; plot - if set, draw chisq vs parameter plot showing 1sigma and 2sigma errors ; ps - if set, create PostScript file of plot ; png - if set, create png file of plot ; nomark - if set, don't draw 1 and 2 sigma lines on plot ; chidelt - don't consider points whose full chisq is > chimin + chidelt (default is 10.) ; connect - if set, connect points with lines (default = 1 if # points to plot is < 50) ; mcfile - name of monte carlo save file to overplot results on chi2 map plot ; ignore_bad - if set, compute and plot sigmas even if bad points were found ; _extra - any keywords to pass to plot command ; ; OUTPUT KEYWORDS: ; out_struct - structure containing results. Tags are: ; name - Parameter number and description:name ; interval - fit interval number ; fit_time - start/end time of fit interval ; param_index - index of parameter chi2 map was done for ; value - value of parameter at minimum chisqr ; s1 - 1-sigma values (2 elements - low side of value, high side of value) ; s2 - 2-sigma values (2 elements - low side of value, high side of value) ; s1_valid - 2 element array of 0s or 1s corresponding to s1. 1 means valid. ; s2_valid - 2 element array of 0s or 1s corresponding to s2. 1 means valid. ; ; Example: ; spex_chi2_map_results, out_struct=out, /plot ; spex_chi2_map_results, /png, savefile='myfile.sav' ; ; OUTPUT: ; The 1- and 2-sigma values for the parameter are printed on the screen. If plot is selected, ; it is drawn to scrren or written in png or ps file. If out_struct passed out, returns ; a summary of the results. ; ; WRITTEN: Kim Tolbert, 23-Aug-2010 ; ; MODIFICATIONS: ; 30-Aug-2010, Kim. Interpolate to find parameter value at chimin+1 and +4. Also don't call ; robust_poly_fit if fewer than 6 points (just find min of chi array in that case). ; 1-Sep-2010, Kim. Changed chimax to chidelt. Added check that have to have at least have 2 trials ; to proceed with finding min and sigmas. Added check whether there are outliers in the derivative - if ; so, can't compute sigmas. ; 19-Oct-2010, Kim. If fewer than 3 points with chi < chimin+chidelt, then set bad, but don't return - ; still do plot. In plot, if bad, set q to finite array, so see all finite points including bad points. ; Added connect keyword; if not set at all, then automatically connect if fewer than 50 points. ; 22-Mar-2011, Kim. Added ignore_bad keyword. ; 05-Apr-2011, Kim. Changed format for output from f9.3 to g12.4 ; 22-Jul-2011, Kim. ; 1. savefile keyword can now contain an array of .sav filenames, and we now loop through ; them all. If writing PS file, all plots written into one save file. ; 2. The reverse histogram method of finding the min doesn't work when there's not a clear, single min. Also ; hsi_outlier doesn't work for data with a slope. Removed those two tests. Now find lowest 5% of points ; (but at least 7 points) and use robust_polyfit on them to find min. ; 3. Added mcfile keyword. If present, and points to a valid ospex monte carlo analysis save file, ; (and one of plot options is selected) the mode, 67% and 95% confidence values from the monte carlo ; analysis will be overplotted on the chi2 mapping plot for the corresponding parameter. ; 8-Oct-2011, Kim. Call al_legend instead of legend (IDL 8. name conflict) ; ;- pro spex_chi2_map_results, savefile=savefile, plot=plot, ps=ps, png=png, nomark=nomark, chidelt=chidelt, $ connect=connect, mcfile=mcfile, out_struct=out_struct, ignore_bad=ignore_bad, _extra=_extra checkvar, chidelt, 10. checkvar, nomark,0 checkvar, outfile, file_basename(savefile) !p.charsize = 1. plot = keyword_set(plot) ps = keyword_set(ps) png = keyword_set(png) if ps then plot=1 if png then plot=1 checkvar, savefile, 'chi2_map.sav' ; show_mc will be 1 if mcfile exists and points to a valid ospex monte carlo analysis save file if keyword_set(mcfile) then begin spex_monte_carlo_results,savefile=mcfile,out_struct=mc if is_struct(mc) then show_mc = 1 endif else show_mc = 0 ; Loop though input chi2_map save files for ifile = 0,n_elements(savefile)-1 do begin infile = savefile[ifile] if ~ file_test(infile, /read) then begin message, 'Can not find or read save file ' + infile, /cont return endif restore, infile bad = n_elements(savechi) lt 2 comp = str2arr(fit_function,'+') for i=0,n_elements(comp)-1 do begin fit_struct = fit_model_components(comp[i], /struct) param_desc = append_arr(param_desc, comp[i] + ': ' + fit_struct.param_desc) endfor nparam = n_elements(savestart) ; # of parameters nf = average(savefullchi/savechi) ; # deg. of freedom out_struct = -1 ind = param_index name = 'Parameter ' + trim(ind) + ': ' + param_desc[ind] finite = where(finite(savefullchi),nfinite) if nfinite eq 0 then begin message, /cont, 'There are no finite values of chisq. Aborting' return endif q = where (savefullchi lt min(savefullchi[finite])+chidelt, nq) if nq lt 3 then begin message, /cont, '< 3 points with chisq < min + chidelt (' + trim(min(savefullchi[finite])+chidelt) + ')' bad=1 endif ; 6/22/2011, remove this test hsi_outlier is not intended for data with a slope - it's wrong ; if not already bad, check if there are any outliers in the relative derivative. If so, don't have ; smooth curve and can't find sigmas. ; if ~bad then begin ; y = savefullchi[q] ; x = savep[ind,q] ; outliers = hsi_outlier(abs(f_div(y(1:*)-y,y)), nsig=3) ; bad = outliers[0] ne -1 ; if bad then message, /cont, 'Data not smooth. Found ' + trim(n_elements(outliers)) + ' outlier(s) in derivative.' ; endif title = 'Interval ' + trim(interval) + ' ' + format_intervals(fit_time,/ut) if keyword_set(ignore_bad) then bad = 0 if ~bad then begin ; find min chisq excluding outliers. First histogram with reverse indices ; to find indices where most points lie (will be at bottom of curve), then ; use robust_poly_fit to fit to that curve excluding outliers, and find ; minimum of the fitted curve. ; h=histogram(savefullchi[q],bins=.2,location=xh,reverse=r) ; dum = max(h,maxh) ; w = r[r[maxh]:r[maxh+1]-1] ; if n_elements(w) gt 6 then begin ; c=robust_poly_fit(savep[ind,q[w]],savefullchi[q[w]],2,yfit) ; dum = min(yfit, miny) ; minind = q[w[miny]] ; endif else begin ; dum = min(savefullchi[q], miny) ; minind = q[miny] ; endelse ; Removed above test. If no single or clear minimum then finding where the most points lie ; doesn't necessarily give you the bottom of the curve. Instead, find the lowest 5% of the values and use ; them to do robust_poly_fit on. s = sort(savefullchi[q]) nuse = .05 * nq > 7 < nq s = q[s[0:nuse-1]] if nuse gt 6 then begin ; robust_poly_fit doesn't work well unless have > 6 points c = robust_poly_fit(savep[ind,s], savefullchi[s], 2, yfit, sig) print,'sig from robust_poly_fit = ' , sig dum = min(yfit,miny) minind = s[miny] endif else begin dum = min(savefullchi[s],miny) minind = s[miny] endelse ;minind is the index of the minimum chisq value, xmin is parameter value there. chimin = savefullchi[minind] xmin = savep[ind,minind] x1 = spex_chi2_map_calcs(savep[ind,*], savefullchi, minind, 1., valid=s1_valid) x2 = spex_chi2_map_calcs(savep[ind,*], savefullchi, minind, 4., valid=s2_valid) sig = !d.name eq 'PS' ? '!9s!3' : '!4r!3' del = !d.name eq 'PS' ? '!9D!3' : '!4D!3' f = '(g12.4)' s1 = [-1.*(xmin-x1[0]), x1[1]-xmin] as1 = arr2str(trim(s1,f),', ') s2 = [-1.*(xmin-x2[0]), x2[1]-xmin] ; 2-sigma values use min chisqr plus 4 as2 = arr2str(trim(s2,f),', ') leg = ['Fit function: ' + fit_function, $ 'Min chisq at ' + trim(xmin,f), $ '1 sigma ' + as1, $ '2 sigma ' + as2] out_struct = {name:name, interval: interval, fit_time: fit_time, param_index: param_index, $ value: xmin, s1: s1, s2: s2, s1_valid: s1_valid, s2_valid: s2_valid} endif else leg = ['Fit function: ' + fit_function, 'Could not determine sigmas!!!'] out = [' ',title, name, leg] prstr,out,/nomore if plot then begin bw = 255 linecolors if ifile eq 0 then dev_save = !d.name if ps and ifile eq 0 then begin ; construct PS output file name. Same as infile, but with .ps extension psfile = file_basename(infile,'.sav',/fold) + '.ps' ps, psfile, /color, /land !p.thick=5 !p.charsize=1. !p.charthick=3 !p.font=0 bw = 0 endif if png then begin pngfile = file_basename(infile,'.sav',/fold) + '.png' tvlct, /get, rr,gg,bb set_plot,'Z' device, set_resolution=[800,500] !p.charsize = .9 endif if ~ps and ~png then window,/free ytitle = 'Full Chisquare' if bad then q = finite ; plot all finite points if there was no clear minimum (bad=1) connect = exist(connect) ? connect : n_elements(q) lt 50 ; connect points if requested, or n<50 psym = connect ? -2 : 2 plot, savep[ind,q],savefullchi[q], psym=psym,symsiz=.3,/yno,$ title=title, ytitle=ytitle, xtitle=name, _extra=_extra if not keyword_set(nomark) then begin if ~bad then begin oplot, [xmin,xmin], [0,1000], color=2 ; vertical line at min value of chisq oplot, !x.crange, [chimin+1, chimin+1],lines=2, color=7 ; horizonatl and vertical lines at chisq min+1 if x1[0] ne 0. then oplot, [x1[0],x1[0]],!y.crange,color=7 if x1[1] ne 0. then oplot, [x1[1],x1[1]],!y.crange,color=7 oplot, !x.crange, [chimin+4, chimin+4],lines=2, color=9 ; horizonatl and vertical lines at chisq min+4 if x2[0] ne 0. then oplot, [x2[0],x2[0]],!y.crange,color=9 if x2[1] ne 0. then oplot, [x2[1],x2[1]],!y.crange,color=9 color = [bw, 2, 7, 9] endif ; if showing monte_carlo results, oplot them in dark blue, solid for mode, dashed for CIs if keyword_set(show_mc) then begin q = where (mc.name eq param_desc[ind], c) if c gt 0 then begin q = q[0] oplot, fltarr(2)+mc[q].mode, !y.crange, color=10 oplot, fltarr(2)+mc[q].v67[0], !y.crange, color=10, linestyle=2 oplot, fltarr(2)+mc[q].v67[1], !y.crange, color=10, linestyle=2 oplot, fltarr(2)+mc[q].v95[0], !y.crange, color=10, linestyle=2 oplot, fltarr(2)+mc[q].v95[1], !y.crange, color=10, linestyle=2 leg = [leg, 'Monte Carlo mode, 67% and 95% CI'] color = [color, 10] endif endif al_legend, leg, textcolor=color, box=0, /center, /top timestamp, /bottom endif if png then begin image = tvrd() write_png, pngfile, image, rr,gg,bb print,'Plot written in file ' + pngfile endif if ifile eq n_elements(savefile)-1 then begin if ps then psclose cleanplot, /silent set_plot, dev_save endif endif ; end of if plot statement endfor ; end of for loop over files message,/cont,'NOTE: Do not trust sigmas without looking at plots (call spex_chi2_map_results with /plot).' end