forked from jdeast/EXOFASTv2
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathexofast_plotchains.pro
167 lines (139 loc) · 6.15 KB
/
exofast_plotchains.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
pro plotchain, pars, label=label, unit=unit, latex=latex, logname=logname, burnndx=burnndx
nallsteps = n_elements(pars[*,0])
nchains = n_elements(pars[0,*])
trimpars = pars[burnndx:*,*]
nsteps = n_elements(trimpars)
medndx = (nsteps/2d0)
halfsigma = erf(1d/sqrt(2d0))/2d
lowsigndx = round( nsteps/2.d0 - nsteps*halfsigma )
hisigndx = round( nsteps/2.d0 + nsteps*halfsigma )
charsize = 1.5
;; check for bad values
bad = where(~finite(trimpars),complement=good)
if bad[0] ne -1 then begin
printandlog, "ERROR: NaNs in " + label + " distribution", logname
return
endif
;; if angular, center distribution about the mode
if unit eq 'DEGREES' then halfrange = 180d0 $
else if unit eq 'RADIANS' then halfrange = !dpi
if unit eq 'DEGREES' or unit eq 'RADIANS' then begin
;; find the mode
hist = histogram(trimpars,nbins=100,locations=x,/nan)
max = max(hist,modendx)
mode = x[modendx]
toohigh = where(pars gt (mode + halfrange))
if toohigh[0] ne -1 then pars[toohigh] -= 2.d0*halfrange
toolow = where(pars lt (mode - halfrange))
if toolow[0] ne -1 then pars[toolow] += 2.d0*halfrange
endif
;; 68% confidence interval
sorted = sort(trimpars)
medvalue = trimpars[sorted[medndx]]
upper = trimpars[sorted[hisigndx]] - medvalue
lower = medvalue - trimpars[sorted[lowsigndx]]
xmax = (medvalue + 4*upper) < max(pars)
xmin = (medvalue - 4*lower) > min(pars)
if xmin eq xmax then begin
printandlog, 'WARNING: ' + label + ' is singularly valued.',logname
endif else begin
;; plot labels
xtitle='!3 Chain link number'
ytitle='!3' + exofast_textoidl(latex)
plot, [0],[0], psym=10, xtitle=xtitle, ytitle=ytitle,$
charsize=charsize,xstyle=1, yrange=[xmin,xmax], xrange=[0,nallsteps],font=1
;; thin the chain for plotting
;; (otherwise the file becomes unusably large)
nplot = 1d3 < n_elements(pars[*,0])
thin = n_elements(pars[*,0])*lindgen(nplot)/(nplot)
for l=0L, nchains-1L do $
oplot, thin, pars[thin,l], color=l*255d0/nchains
oplot, [burnndx,burnndx],[-9d99,9d99]
endelse
end
pro exofast_plotchains, ss, chainfile=chainfile, logname=logname
if n_elements(chainfile) eq 0 then chainfile = 'chain.ps'
chi2 = reform((*ss.chi2),ss.nsteps/ss.nchains,ss.nchains)
burnndx = ss.burnndx
goodchains = (*ss.goodchains)
;burnndx = getburnndx(chi2,goodchains=goodchains, logname=logname)
;; prepare the postscript device
mydevice=!d.name
set_plot, 'PS'
aspect_ratio=1.5
xsize = 18
ysize=xsize/aspect_ratio
!p.font=0
defsysv, '!GDL', exists=runninggdl
if runninggdl then chainfile0 = file_dirname(chainfile) + path_sep() + file_basename(chainfile,'.ps') + '.001.ps' $
else chainfile0 = chainfile
device, filename=chainfile0, set_font='Times',/tt_font, xsize=xsize,ysize=ysize,/color,bits=24
loadct,39,/silent
red = 254
!p.multi=[0,2,4] ;; 8 to a page
loglike = -chi2/2d0
ymin = min(loglike,max=ymax)
if ymin lt 0 then sign = ' + ' else sign = ' - '
latex = 'log(Like)' + sign + strtrim(abs(ymax),2)
plotchain, loglike[*,goodchains]-ymax, latex=latex, unit='', label='loglike', logname=logname, burnndx=burnndx
page = 2
for i=0, n_tags(ss)-1 do begin
for j=0, n_elements(ss.(i))-1 do begin
for k=0, n_tags(ss.(i)[j])-1 do begin
derive = 0B
m=-1
;; this captures the detrending variables
if (size(ss.(i)[j].(k)))[1] eq 10 then begin
if ptr_valid(ss.(i)[j].(k)) then begin
for l=0L, n_tags(*(ss.(i)[j].(k)))-1 do begin
if (size((*(ss.(i)[j].(k))).(l)))[2] eq 8 then begin
for m=0L, n_elements((*(ss.(i)[j].(k))).(l))-1 do begin
if tag_exist((*(ss.(i)[j].(k))).(l)[m],'derive') then begin
if (*(ss.(i)[j].(k))).(l)[m].derive or (*(ss.(i)[j].(k))).(l)[m].fit then begin
;; GDL can't do multi-page plots
if runninggdl and !p.multi[0] eq 0 then begin
device, /close
chainfile0 = file_dirname(chainfile) + path_sep() + file_basename(chainfile,'.ps') + '.' + string(page,format='(i03)') + '.ps'
device, filename=chainfile0, set_font='Times',/tt_font, xsize=xsize,ysize=ysize,/color,bits=24
page += 1
endif
pars = (reform((*(ss.(i)[j].(k))).(l)[m].value,ss.nsteps/ss.nchains,ss.nchains))[*,goodchains]
plotchain, pars, label=(*(ss.(i)[j].(k))).(l)[m].label,$
unit = (*(ss.(i)[j].(k))).(l)[m].unit,$
latex = (*(ss.(i)[j].(k))).(l)[m].latex,$
logname=logname, burnndx=burnndx
endif
endif
endfor
endif
endfor
endif
endif else if n_tags(ss.(i)[j].(k)) ne 0 then begin
if n_tags(ss.(i)[j].(k)) ne 0 then begin
if tag_exist(ss.(i)[j].(k),'derive') then begin
if ss.(i)[j].(k).derive then begin
;; GDL can't do multi-page plots
if runninggdl and !p.multi[0] eq 0 then begin
device, /close
chainfile0 = file_dirname(chainfile) + path_sep() + file_basename(chainfile,'.ps') + '.' + string(page,format='(i03)') + '.ps'
device, filename=chainfile0, set_font='Times',/tt_font, xsize=xsize,ysize=ysize,/color,bits=24
page += 1
endif
pars = (reform(ss.(i)[j].(k).value,ss.nsteps/ss.nchains,ss.nchains))[*,goodchains]
plotchain, pars, label=ss.(i)[j].(k).label,$
unit = ss.(i)[j].(k).unit,$
latex = ss.(i)[j].(k).latex,$
logname=logname, burnndx=burnndx
endif
endif
endif
endif
endfor
endfor
endfor
;; clean up the postscript device
device, /close
!p.font=-1
!p.multi=0
set_plot, mydevice
end