forked from jdeast/EXOFASTv2
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathexofast_errell.pro
executable file
·191 lines (173 loc) · 6.92 KB
/
exofast_errell.pro
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
;+
; NAME:
; exofast_errell
;
; PURPOSE:
; Gives the path of the error ellipses at constant probability given
; two parameter distributions.
;
; PROCEDURE
; 1. Creates a 2D histogram of the data with NBINS bins per
; parameter.
; 2. Calculates the fraction of data points above each value in the
; histogram.
; 3. Compares each calculated fraction with the probabilities,
; and defines the closest as the levels to pass to CONTOUR.
; 4. Creates a contour plot of the 2D histogram
; 5. Scales the path to be in the same units as the input and
; recreates the path, adding NaNs between each level for easy
; plotting
;
; CALLING SEQUENCE:
; exofast_errell, xpar, ypar [,NXBIN=nxbin,NYBIN=nybin,PROBS=PROBS]$
; [,XPATH=xpath][,YPATH=ypath][,/PLOT][,/OPLOT]
;
; INPUTS:
; XPAR - A vector of parameters to be plotted on the x axis.
; YPAR - A vector of parameters to be plotted on the y axis.
;
; OPTIONAL INPUTS:
; NXBIN - The number of bins to take in each direction. Default
; is 100. If you have only a small number of values for
; xpar and ypar (10000), this should be smaller.
; NYBIN - The number of bins to take in each direction. Default
; is NXBIN.
; XMIN/YMIN - The minimum x/y value to make the contours
; XMAX/YMAX - The maximum x/y value to make the contours
; PROBS - An array containing the probability contours to
; plot. The default is
; probs = erf([1.d0,2.d0]/sqrt(2.d0)), or
; probs = [0.6826894921370859d0,$
; 0.9544997361036416d0]
; Note: to have a meaningful contour enclosing 99.7% of
; the points, you must have ~1 million points or make
; NBINS smaller.
;
; OPTIONAL KEYWORDS:
; PLOT - Plot the contours on a new plot
; OPLOT - Overplot the contours on an existing plot
;
; OPTIONAL OUTPUTS:
; XPATH/YPATH - The x/y values of the contours. For more
; customizable plotting, return these and plot inside
; your own program. All contours can be plotted simply
; with IDL> plot, xpath, ypath
; NOTE: XPATH/YPATH contain NaNs between levels to make the
; plotting easier.
; OUTSIDENDX - The indices of the parameters outside of the
; outer-most contour (requires David Fanning's Inside
; procedure)
;
; EXAMPLE:
; plot the 68% and 95% contours, and every point outside them
;
; x = randomn(seed,100000)
; y = randomn(seed,100000)
; exofast_errell, x,y,out=out,probs=[0.68,0.95],xpath=xpath,ypath=ypath
; plot, xpath, ypath, /iso
; if out[0] ne -1 then oplot, x[out],y[out],psym=3
;
; DEPENDENCIES:
; INSIDE (Coyote library):
; http://www.idlcoyote.com/programs/inside.pro
;
; NOTE: there is non-vectorized version of his code online, which will
; not work. Use the one at the link above.
;
; REVISION HISTORY:
; 2010/03/01 - Written: Jason Eastman - The Ohio State University
; 2010/11/05 - Added OUTSIDENDX output (inspired by Ben Shappee, OSU)
; Added, NXBIN, XMIN, XMAX, NYBIN, YMIN, YMAX inputs.
; 2011/03/25 - Better memory management for INSIDE -- split up
; array to check (takes a little longer, but actually
; works with large arrays!)
;-
pro exofast_errell, xpar, ypar, probs=probs, xpath=xpath, ypath=ypath, $
plot=plot, oplot=oplot, outsidendx=outsidendx, $
nxbin=nxbin, xmin=xmin, xmax=xmax, $
nybin=nybin, ymin=ymin, ymax=ymax,logname=logname
;; 1 and 2 sigma contours (as defined in 1D)
if n_elements(probs) eq 0 then probs = erf([1.d0,2.d0]/sqrt(2.d0))
;; define the bounds of the histogram in the x direction
if n_elements(xmin) eq 0 then xmin = min(xpar,/nan)
if n_elements(xmax) eq 0 then xmax = max(xpar,/nan)
if n_elements(nxbin) eq 0 then nxbin = 100
if xmin eq xmax then begin
printandlog, 'Histogram has zero width (parameter has no dispersion?)',logname
return
endif
xbinsz = double(xmax-xmin)/(nxbin-1)
;; define the bounds of the histogram in the y direction
if n_elements(ymin) eq 0 then ymin = min(ypar,/nan)
if n_elements(ymax) eq 0 then ymax = max(ypar,/nan)
if n_elements(nybin) eq 0 then nybin = nxbin
if ymin eq ymax then begin
printandlog, 'Histogram has zero width (parameter has no dispersion?)',logname
return
endif
ybinsz = double(ymax-ymin)/(nybin-1)
;; create the 2D histogram
hist2d = hist_2d(xpar,ypar,$
bin1=xbinsz,min1=xmin,max1=xmax,$
bin2=ybinsz,min2=ymin,max2=ymax)
;; find the fraction of points enclosed above each bin
;; NOTE: includes points excluded by xmin, xmax, ymin, ymax
maxhist = max(hist2d)
enclosed = dblarr(maxhist+1)
for i=0L, maxhist do begin
good = where(hist2d ge i)
if good[0] ne -1 then enclosed[i] = total(hist2d[good])
endfor
fraction = enclosed/n_elements(where(finite(xpar)))
;; find the levels that correspond to the different values of sigma
levels = dblarr(n_elements(probs))
for i=0, n_elements(probs)-1 do begin
min = min(abs(fraction-probs[i]),level)
levels[i] = level
endfor
levels = reverse(levels)
;; make sure LEVELS is a sensible input to CONTOUR
if n_elements(uniq(levels)) ne n_elements(probs) then begin
printandlog, 'Resolution too fine -- lower NXBIN or use more elements!',logname
return
endif
;; get the paths of the contours
contour, hist2d, xtitle=xtitle,ytitle=ytitle, levels=levels,$
path_info=info,path_xy=xy,/path_data_coords,/overplot,/path_double
;; prep for outside index
if arg_present(outsidendx) then isinside = lonarr(n_elements(xpar))
;; convert the path values input units
x = xmin + (xy[0,*]+1)*(xmax-xmin)/nxbin
y = ymin + (xy[1,*]+1)*(ymax-ymin)/nybin
;; convert the path to intuitively plottable path
;; ie, plot, xpath, ypath
xpath = !values.d_nan
ypath = !values.d_nan
for i=0, n_elements(info)-1 do begin
S = [indgen(info(i).n), 0]
xpath = [xpath,x[info(i).offset + S],!values.d_nan]
ypath = [ypath,y[info(i).offset + S],!values.d_nan]
;; finds elements that are outside all outer-level contours
;; (isinside = 0)
if info(i).level eq 0 and arg_present(outsidendx) then begin
check = where(isinside eq 0, ncheck)
if check[0] ne -1 then begin
;; inside is very memory inefficient --
;; split up the check to more managable bites
for j=0L, ncheck, 1d4 do begin
checknow = j + indgen(1d4 < (ncheck-j))
isinside[check[checknow]] += $
inside(xpar[check[checknow]],ypar[check[checknow]],$
x[info(i).offset+S],y[info(i).offset+S])
endfor
endif
endif
endfor
if arg_present(outsidendx) then outsidendx = where(isinside eq 0)
;; begin a new plot if desired
if keyword_set(plot) then $
plot, xpath, ypath,_extra=extra
;; overplot if desired
if keyword_set(oplot) then $
oplot, xpath, ypath
end