diff --git a/derivepars.pro b/derivepars.pro index 82bcc9d35..e829b2639 100644 --- a/derivepars.pro +++ b/derivepars.pro @@ -61,18 +61,18 @@ for i=0, ss.nplanets-1 do begin ss.planet[i].period.value[positive],$ ss.star.mstar.value[positive]) - ss.planet[i].mp.value = ss.planet[i].mpsun.value/mjup - ss.planet[i].msini.value = ss.planet[i].mp.value*sini - ss.planet[i].mpearth.value = ss.planet[i].mpsun.value/mearth - ss.planet[i].msiniearth.value = ss.planet[i].mpearth.value*sini - ss.planet[i].q.value = ss.planet[i].mpsun.value/ss.star.mstar.value + ss.planet[i].mp.value = ss.planet[i].mpsun.value/mjup ;; m_jupiter + ss.planet[i].msini.value = ss.planet[i].mp.value*sini ;; m_jupiter + ss.planet[i].mpearth.value = ss.planet[i].mpsun.value/mearth ;; m_earth + ss.planet[i].msiniearth.value = ss.planet[i].mpearth.value*sini ;; m_earth + ss.planet[i].q.value = ss.planet[i].mpsun.value/ss.star.mstar.value ;; unitless - ss.planet[i].arsun.value=(G*(ss.star.mstar.value+ss.planet[i].mp.value)*ss.planet[i].period.value^2/(4d0*!dpi^2))^(1d0/3d0) - ss.planet[i].ar.value = ss.planet[i].arsun.value/ss.star.rstar.value ;; unitless - ss.planet[i].a.value = ss.planet[i].arsun.value/AU ;; AU - ss.planet[i].rpsun.value = ss.planet[i].p.value*ss.star.rstar.value - ss.planet[i].rp.value = ss.planet[i].rpsun.value/rjup - ss.planet[i].rpearth.value = ss.planet[i].rpsun.value/rearth + ss.planet[i].arsun.value=(G*(ss.star.mstar.value+ss.planet[i].mpsun.value)*ss.planet[i].period.value^2/(4d0*!dpi^2))^(1d0/3d0) ;; (a1 + a2)/rsun + ss.planet[i].ar.value = ss.planet[i].arsun.value/ss.star.rstar.value ;; (a1 + a2)/rstar + ss.planet[i].a.value = ss.planet[i].arsun.value/AU ;; AU + ss.planet[i].rpsun.value = ss.planet[i].p.value*ss.star.rstar.value ;; r_sun + ss.planet[i].rp.value = ss.planet[i].rpsun.value/rjup ;; r_jupiter + ss.planet[i].rpearth.value = ss.planet[i].rpsun.value/rearth ;; r_earth ;; time of periastron ss.planet[i].phase.value=exofast_getphase(ss.planet[i].e.value,ss.planet[i].omega.value,/pri) @@ -151,12 +151,12 @@ for i=0, ss.nplanets-1 do begin ss.planet[i].ps.value = (ss.star.rstar.value-ss.planet[i].rpsun.value)/ss.planet[i].arsun.value*(1d0 - ss.planet[i].esinw.value)/(1d0-ss.planet[i].e.value^2) ;; eq 9, Winn 2010 ss.planet[i].rhop.value = ss.planet[i].mpsun.value/(ss.planet[i].rpsun.value^3)*1.41135837d0 - ss.planet[i].loggp.value = alog10(G*ss.planet[i].mp.value/ss.planet[i].rp.value^2*9.31686171d0) ;; cgs + ss.planet[i].loggp.value = alog10(G*ss.planet[i].mpsun.value/ss.planet[i].rpsun.value^2*9.31686171d0) ;; cgs ss.planet[i].safronov.value = ss.planet[i].ar.value*ss.planet[i].q.value/ss.planet[i].p.value ;; depth != delta if grazing (ignore limb darkening) ss.planet[i].delta.value = ss.planet[i].p.value^2 - for j=0, ss.nsteps-1 do begin + for j=0L, ss.nsteps-1L do begin exofast_occultquad, ss.planet[i].b.value[j], 0d0, 0d0, ss.planet[i].p.value[j],mu1 ss.planet[i].depth.value[j] = 1d0-mu1[0] endfor diff --git a/exofast_chi2v2.pro b/exofast_chi2v2.pro index 6a20f5512..9144c8c7e 100644 --- a/exofast_chi2v2.pro +++ b/exofast_chi2v2.pro @@ -246,10 +246,12 @@ ss.star.mstar.value = 10^ss.star.logmstar.value ;; use the YY tracks to guide the stellar parameters if ~ss.noyy then begin if keyword_set(psname) then begin - chi2 += massradius_yy3(ss.star.mstar.value, ss.star.feh.value, ss.star.age.value, ss.star.teff.value,yyrstar=rstar, debug=ss.debug, psname=psname+'.yy.eps') + yychi2 = massradius_yy3(ss.star.mstar.value, ss.star.feh.value, ss.star.age.value, ss.star.teff.value,yyrstar=rstar, debug=ss.debug, psname=psname+'.yy.eps') endif else begin - chi2 += massradius_yy3(ss.star.mstar.value, ss.star.feh.value, ss.star.age.value, ss.star.teff.value,yyrstar=rstar, debug=ss.debug) + yychi2 = massradius_yy3(ss.star.mstar.value, ss.star.feh.value, ss.star.age.value, ss.star.teff.value,yyrstar=rstar, debug=ss.debug) endelse + chi2 += yychi2 + ;print, 'YY penalty = ' + strtrim(yychi2,2) if ~finite(chi2) then begin if ss.debug then print, 'star is bad' return, !values.d_infinity @@ -263,6 +265,11 @@ if ss.star.errscale.value le 0 then begin return, !values.d_infinity endif +if ss.star.alpha.value lt -0.3d0 or ss.star.alpha.value gt 0.7d0 then begin + if ss.debug then print, 'alpha is bad (' + strtrim(ss.star.alpha.value,2) + ')' + return, !values.d_infinity +endif + if ss.amoeba then begin ;; need some derived parameters ;; (but others will change after we change logk or p, so @@ -437,11 +444,11 @@ if file_test(ss.star.fluxfile) then begin if keyword_set(psname) then begin sedchi2 = exofast_sed(ss.star.fluxfile, ss.star.teff.value, ss.star.rstar.value,$ ss.star.av.value, ss.star.distance.value, $ - logg=ss.star.logg.value,met=ss.star.feh.value,verbose=ss.debug, f0=f, fp0=fp, ep0=ep, psname=psname+'.sed.eps') + logg=ss.star.logg.value,met=ss.star.feh.value,alpha=ss.star.alpha.value,verbose=ss.debug, f0=f, fp0=fp, ep0=ep, psname=psname+'.sed.eps') endif else begin sedchi2 = exofast_sed(ss.star.fluxfile, ss.star.teff.value, ss.star.rstar.value,$ ss.star.av.value, ss.star.distance.value, $ - logg=ss.star.logg.value,met=ss.star.feh.value,verbose=ss.debug, f0=f, fp0=fp, ep0=ep) + logg=ss.star.logg.value,met=ss.star.feh.value,alpha=ss.star.alpha.value,verbose=ss.debug, f0=f, fp0=fp, ep0=ep) endelse if ~finite(sedchi2) then begin @@ -455,6 +462,7 @@ if file_test(ss.star.fluxfile) then begin return, !values.d_infinity endif chi2 += sedchi2 +; print, 'SED penalty = ' + strtrim(sedchi2,2) endif ;; RV model (non-interacting planets) @@ -502,6 +510,7 @@ for j=0, ntelescopes-1 do begin if ~finite(rvchi2) then stop chi2 += rvchi2 +; print, 'RV penalty = ' + strtrim(rvchi2,2) ; chi2 += total(((rv.rv - modelrv)/rv.err)^2) endfor @@ -532,6 +541,8 @@ for j=0, ntransits-1 do begin u2err = 0.05d0 chi2 += ((band.u1.value-u1claret)/u1err)^2 chi2 += ((band.u2.value-u2claret)/u2err)^2 +; print, 'u1 penalty = ' + strtrim(((band.u1.value-u1claret)/u1err)^2,2) +; print, 'u2 penalty = ' + strtrim(((band.u2.value-u2claret)/u2err)^2,2) endif ;; Kepler Long candence data; create several model points and average @@ -547,6 +558,14 @@ for j=0, ntransits-1 do begin modelflux = dblarr(npoints) + 1d0 endelse + ;; get the motion of the star due to the planet + junk = exofast_getb2(transitbjd,inc=ss.planet.i.value,a=ss.planet.ar.value,$ + tperiastron=ss.planet.tp.value,$ + period=ss.planet.period.value,$ + e=ss.planet.e.value,omega=ss.planet.omega.value,$ + q=ss.star.mstar.value/ss.planet.mpsun.value,$ + x1=x1,y1=y1,z1=z1) + for i=0, nplanets-1 do begin if ss.planet[i].fittran then begin @@ -568,8 +587,8 @@ for j=0, ntransits-1 do begin reflect=band.reflect.value, $ dilute=band.dilute.value,$ tc=ss.planet[i].tc.value,$ - rstar=ss.star.rstar.value/AU) - 1d0) - + rstar=ss.star.rstar.value/AU,$ + x1=x1,y1=y1,z1=z1) - 1d0) endif endfor @@ -584,6 +603,7 @@ for j=0, ntransits-1 do begin ;; chi^2 transitchi2 = exofast_like(transit.flux - modelflux,ss.transit[j].variance.value,transit.err,/chi2) + (*ss.transit[j].transitptrs).residuals = transit.flux - modelflux (*ss.transit[j].transitptrs).model = modelflux @@ -593,11 +613,14 @@ for j=0, ntransits-1 do begin forprint, transitbjd, modelflux, format='(f0.8,x,f0.6)', textout=psname + '.model.' + strtrim(j,2) + '.flux.txt', /nocomment,/silent if ~finite(transitchi2) then stop - chi2 += transitchi2 -; print, 'transit penalty: ' + strtrim(transitchi2,2) -; chi2 += total(((transit.flux - modelflux)/transit.err)^2) +;transitchi2 = strtrim(total(((transit.flux - modelflux)/transit.err)^2),2) + + chi2 += transitchi2 +; print, 'transit penalty: ' + strtrim(transitchi2,2) + ' ' + strtrim(ss.transit[j].variance.value,2) +;print, 'transit chi2 = ' + strtrim(total(((transit.flux - modelflux)/transit.err)^2),2) +;stop ; screen = GET_SCREEN_SIZE() ; if win_state(5) then wset, 5 $ ; else window, 5, xsize=xsize, ysize=ysize, xpos=screen[0]/3d0, ypos=0 @@ -638,6 +661,9 @@ if ss.ttvs then begin chi2 += ((coeffs[0]-ss.planet[i].tc.value)/sigma[0])^2 chi2 += ((coeffs[1]-ss.planet[i].period.value)/sigma[1])^2 +; print, ((coeffs[0]-ss.planet[i].tc.value)/sigma[0])^2 +; print, ((coeffs[0]-ss.planet[i].tc.value)/sigma[0])^2 + if ss.debug or keyword_set(psname) then begin if keyword_set(psname) then begin ;; astrobetter.com tip on making pretty IDL plots @@ -727,7 +753,6 @@ if ss.debug then print, ss.star.rstar.value, pars, chi2, format='(' + strtrim(n_ ;; if this stop is triggered, you've found a bug!! if ~finite(chi2) then stop - return, chi2 end diff --git a/exofast_demc.pro b/exofast_demc.pro index 24797aa92..7e09531e6 100644 --- a/exofast_demc.pro +++ b/exofast_demc.pro @@ -155,7 +155,7 @@ pro exofast_demc, bestpars,chi2func,pars,chi2=chi2, tofit=tofit,$ nchains=nchains, angular=angular,$ burnndx=burnndx, removeburn=removeburn, $ gelmanrubin=gelmanrubin, tz=tz, maxgr=maxgr, mintz=mintz, $ - derived=derived, resumendx=resumendx + derived=derived, resumendx=resumendx,stretch=stretch ;; default values if n_elements(maxsteps) eq 0 then maxsteps = 1d5 @@ -223,8 +223,9 @@ for j=0, nchains-1 do begin ;; start 5 steps from best value. If that's not allowed ;; (infinite chi^2), slowly approach 0 steps away from the best value ;; niter should never be larger than ~5000 + if j eq 0 then factor = 0d0 else factor = 1d0 ;; first chain starts at best fit pars[0:nfit-1,0,j] = bestpars[tofit] + $ - 5d0/exp(niter/1000d)*scale*call_function(randomfunc,seed,nfit,/normal) + 5d0/exp(niter/1000d)*scale*call_function(randomfunc,seed,nfit,/normal)*factor newpars[tofit] = pars[0:nfit-1,0,j] ;; find the chi^2 chi2[0,j] = call_function(chi2func, newpars, determinant=det, derived=dpar) @@ -247,6 +248,8 @@ print if n_elements(dpar) ne 0 then oldderived = reform(derived[*,0,*],nderived,nchains) +a = 2d0 ;; a > 1 + nextrecalc = 100L npass = 1L nstop = 0L @@ -255,6 +258,9 @@ tz0 = 0d0 tzsteps = 0L alreadywarned = 0L t0 = systime(/seconds) +minburnndx = 0L + +idealacceptancerate = 0.23d0 ;; automatically scale "a" to get ideal acceptance rate ;; start MCMC chain for i=resumendx,maxsteps-1L do begin @@ -266,18 +272,38 @@ for i=resumendx,maxsteps-1L do begin ;; automatically thin the chain (saves memory) for k=0L, nthin-1L do begin - ;; differential evolution mcmc step -- see Ter Braak 2006 - repeat r1=floor(call_function(randomfunc,seed)*(nchains)) $ - until r1 ne j ;; a random chain - repeat r2=floor(call_function(randomfunc,seed)*(nchains)) $ - until r2 ne j and r2 ne r1 ;; another random chain - newpars[tofit] = oldpars[tofit] + $ - gamma*(pars[*,i-1,r1]-pars[*,i-1,r2]) + $ - (call_function(randomfunc,seed,nfit)-0.5d0)*scale/10d0 - + if keyword_set(stretch) then begin + ;; affine invariant "stretch move" step + ;; see Goodman & Weare, 2010 (G10), eq 7 + ;; http://msp.org/camcos/2010/5-1/camcos-v5-n1-p04-p.pdf + repeat r1=floor(call_function(randomfunc,seed)*(nchains)) $ + until r1 ne j ;; a random chain + + ;; a random number between 1/a and a, weighted toward lower + ;; values to satisfy the symmetry condition (G10, eq 8) + z = ((a-1d0)*call_function(randomfunc,seed) + 1d0)^2d0/2d0 + + ;; use the most recent step (Goodman, pg 69) + if r1 lt j then ndx = i else ndx = i-1 + ;; define the new step (eq 7 in G10) + newpars[tofit] = pars[*,ndx,r1] + z*(oldpars[tofit]-pars[*,ndx,r1]) + fac = z^(nfit-1) ;; GW10, last equation, pg 70 + endif else begin + ;; differential evolution mcmc step (Ter Braak, 2006) + repeat r1=floor(call_function(randomfunc,seed)*(nchains)) $ + until r1 ne j ;; a random chain + repeat r2=floor(call_function(randomfunc,seed)*(nchains)) $ + until r2 ne j and r2 ne r1 ;; another random chain + newpars[tofit] = oldpars[tofit] + $ + gamma*(pars[*,i-1,r1]-pars[*,i-1,r2]) + $ + (call_function(randomfunc,seed,nfit)-0.5d0)*scale/10d0 + fac = 1d0 + endelse + + ;; calculate the chi^2 of the new step newchi2 = call_function(chi2func,newpars,determinant=det,derived=dpar) - C = olddet[j]/det*exp((oldchi2 - newchi2)/2d0) + C = fac*olddet[j]/det*exp((oldchi2 - newchi2)/2d0) ;; accept the step; update values if call_function(randomfunc,seed) lt C then begin @@ -378,6 +404,26 @@ for i=resumendx,maxsteps-1L do begin ;; print the progress, cumulative acceptance rate, and time remaining acceptancerate = strtrim(string(naccept/double(i*nchains*nthin)*100,$ format='(f6.2)'),2) + + ;; Limited testing suggests this is unnecessary (counter-productive, even) +; ;; tune the stretch scale (a) to achieve the optimal acceptance rate (~23%) +; ;; if we've taken fewer than 5% of the total steps +; if double(i)/maxsteps lt 0.05 and keyword_set(stretch) then begin +; if naccept/(double(i*nchains*nthin)) lt idealacceptancerate*0.75 or $ +; naccept/(double(i*nchains*nthin)) gt idealacceptancerate*1.25 then begin +; +; print +; print, 'Acceptance rate (' + acceptancerate + ') not ideal -- scaling stretch scale from ' + strtrim(a,2) + ' to ' + strtrim(a*naccept/(double(i*nchains*nthin))/idealacceptancerate,2) +; a = a*naccept/(double(i*nchains*nthin))/idealacceptancerate +; +; ;; more continous scaling +; if accepted then a += (1d0-idealacceptancerate)/i $ +; else a -= idealacceptancerate/i +; +; minburnndx = i +; endif +; endif + timeleft = (systime(/seconds)-t0)*(maxsteps/(i+1d)-1d) units = ' seconds ' if timeleft gt 60 then begin @@ -418,7 +464,7 @@ medchi2 = min(median(chi2[0.1*nstop:nstop,*],dimension=1)) ;; good chains must have some values below this median goodchains = where(minchi2 lt medchi2,ngood) -burnndx = 0L +burnndx = minburnndx for j=0L, ngood-1 do begin tmpndx = (where(chi2[0:nstop,goodchains[j]] lt medchi2))(0) if tmpndx gt burnndx then burnndx = tmpndx diff --git a/exofast_getb.pro b/exofast_getb.pro index ec04856fc..9c354dfdb 100755 --- a/exofast_getb.pro +++ b/exofast_getb.pro @@ -53,9 +53,7 @@ function exofast_getb, bjd, i=i, a=a, tperiastron=tperiastron, Period=P, $ e=e, omega=omega, x=x, y=y, z=z, $ - lonascnode=lonascnode, q=q - -if n_elements(q) eq 0 then q = !values.d_infinity + lonascnode=lonascnode ;; calculate the mean anomaly corresponding to each observed time meananom = (2.d0*!dpi*(1.d0 + (bjd - Tperiastron)/P)) mod (2.d0*!dpi) @@ -71,10 +69,8 @@ endif else begin e=0.d0 endelse -atot = a ;; a1 + a2 - ;; calculate the corresponding (x,y) coordinates of planet -r = atot*(1d0-e^2)/(1d0+e*cos(trueanom)) +r = a*(1d0-e^2)/(1d0+e*cos(trueanom)) ;; as seen from observer x = -r*cos(trueanom + omega) diff --git a/exofast_getmcmcscale.pro b/exofast_getmcmcscale.pro index 26490925c..5264e681d 100755 --- a/exofast_getmcmcscale.pro +++ b/exofast_getmcmcscale.pro @@ -144,7 +144,7 @@ for i=0, nfit-1 do begin ' (' + strtrim(chi2,2) + ')' endif - if (origbestchi2 - bestchi2) gt 0.1 then betterfound = 1B + if (origbestchi2 - bestchi2) gt 1d0 then betterfound = 1B ;; chi2 is actually lower! (Didn't find the best fit) ;; attempt to fix @@ -153,8 +153,9 @@ for i=0, nfit-1 do begin ;; could be way off, double the step for faster convergence mcmcscale[i,j] *= 2d0 bestchi2 = chi2 - chi2changed = 1 niter = 0 + chi2changed = 1 + endelse deltachi2 = chi2-bestchi2 @@ -171,14 +172,17 @@ for i=0, nfit-1 do begin if not chi2changed then begin if abs(bestdeltachi2 - 1.d0) lt 0.75 then begin mcmcscale[i,j] = bestscale - chi2 = bestchi2 + 1 + chi2 = bestchi2 + 1d0 endif else begin - message, 'Convergence Error: cannot find the value' + $ - ' for which deltachi^2 = 1 for parameter ' + strtrim(tofit[i],2),/continue + message, 'Cannot find the value for which deltachi^2 = 1 for parameter ' +$ + strtrim(tofit[i],2) + '; assuming rough chi^2 surface. Using delta chi^2 = ' +$ + strtrim(bestdeltachi2,2) + ' and scaling step (' + strtrim(mcmcscale[i,j],2) +$ + ') to ' + strtrim(mcmcscale[i,j]*(bestchi2+1)/chi2,2),/continue + + ;; extrapolate scale to delta chi^2 = 1 + ;; (may cause issues near boundaries) mcmcscale[i,j] = mcmcscale[i,j]*(bestchi2+1)/chi2 - chi2 = bestchi2 + 1 -;stop - if tofit[i] eq 4 then stop + chi2 = bestchi2 + 1d0 endelse endif endif @@ -199,10 +203,10 @@ for i=0, nfit-1 do begin ;; near a boundary ;; exclude this direction from the scale calculation - if ~finite(chi2) then begin - chi2 = bestchi2 + 1.d0 - mcmcscale[i,j] = minstep - endif +; if ~finite(chi2) then begin +; chi2 = bestchi2 + 1.d0 +; mcmcscale[i,j] = minstep +; endif niter++ ;; print the steps and chi^2 it's getting @@ -215,10 +219,6 @@ for i=0, nfit-1 do begin endif next: -;if tofit[i] eq 4 and j eq 1 and n_elements(tofit) gt 4 then begin -; print, mcmcscale -; stop -;endif endrep until abs(chi2 - bestchi2 - 1.d0) lt 1d-8 endfor endfor @@ -235,11 +235,11 @@ if betterfound and ~keyword_set(skipiter) then begin endif ;; replace undefined errors with the other -bad = where(~finite(mcmcscale[*,0])) +bad = where(~finite(mcmcscale[*,0]) or mcmcscale[*,0] eq 0d0) if bad[0] ne -1 then mcmcscale[bad,0] = mcmcscale[bad,1] -bad = where(~finite(mcmcscale[*,1])) +bad = where(~finite(mcmcscale[*,1]) or mcmcscale[*,1] eq 0d0) if bad[0] ne -1 then mcmcscale[bad,1] = mcmcscale[bad,0] -bad = where(~finite(mcmcscale),nbad) +bad = where(~finite(mcmcscale) or mcmcscale eq 0d0,nbad) if bad[0] ne -1 then return, -1 ;; return the average of high and low diff --git a/exofast_plotdist2.pro b/exofast_plotdist2.pro index 91d99304b..f73db6a42 100755 --- a/exofast_plotdist2.pro +++ b/exofast_plotdist2.pro @@ -65,11 +65,11 @@ if n_elements(covarname) eq 0 then covarname = 'covar.ps' if n_elements(probs) eq 0 then probs = erf([1d,2d]/sqrt(2d0)) burnndx = ss.burnndx -nsteps = ss.nsteps - burnndx -medndx = nsteps/2 +nsteps = ss.nsteps/ss.nchains +medndx = ((nsteps-ss.burnndx)/2d0) * ss.nchains halfsigma = erf(1d/sqrt(2d0))/2d -lowsigndx = round(nsteps/2.d0 - nsteps*halfsigma) -hisigndx = round(nsteps/2.d0 + nsteps*halfsigma) +lowsigndx = round(((nsteps-ss.burnndx)/2.d0 - (nsteps-ss.burnndx)*halfsigma)*ss.nchains) +hisigndx = round(((nsteps-ss.burnndx)/2.d0 + (nsteps-ss.burnndx)*halfsigma)*ss.nchains) ;; prepare the postscript device mydevice=!d.name @@ -80,8 +80,10 @@ ysize=xsize/aspect_ratio !p.font=0 device, filename=pdfname device, set_font='Times',/tt_font -device, xsize=xsize,ysize=ysize +device, xsize=xsize,ysize=ysize,/color,bits=24 charsize=1.5 +loadct,39,/silent +red = 254 allpars = [] parnames = [] @@ -102,7 +104,8 @@ for i=0, n_tags(ss)-1 do 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 then begin - pars = (*(ss.(i)[j].(k))).(l)[m].value[burnndx:*,*] + pars = (reform((*(ss.(i)[j].(k))).(l)[m].value,nsteps,ss.nchains))[burnndx:*,*] +; pars = (*(ss.(i)[j].(k))).(l)[m].value[burnndx:*,*] derive = 1B label = (*(ss.(i)[j].(k))).(l)[m].label unit = (*(ss.(i)[j].(k))).(l)[m].unit @@ -118,7 +121,7 @@ for i=0, n_tags(ss)-1 do 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 - pars = ss.(i)[j].(k).value[burnndx:*,*] + pars = (reform(ss.(i)[j].(k).value,nsteps,ss.nchains))[burnndx:*,*] derive = 1B label = ss.(i)[j].(k).label unit = ss.(i)[j].(k).unit @@ -176,10 +179,22 @@ for i=0, n_tags(ss)-1 do begin ;; plot labels xtitle='!3' + textoidl(latex + '_{' + ss.(i)[j].label + '}') ytitle='!3Probability' - + + + hist = histogram(pars,nbins=100,locations=x,min=xmin,max=xmax) - plot, x, hist/double(nsteps), psym=10, xtitle=xtitle, ytitle=ytitle,$ + plot, x, hist/double(total(hist)), psym=10, xtitle=xtitle, ytitle=ytitle,$ charsize=charsize,xstyle=1, xrange=[xmin,xmax], font=1 + + for l=0L, ss.nchains-1L do begin + hist = histogram(pars[*,l],nbins=100,locations=x,min=xmin,max=xmax) + if total(hist) gt 0 then $ + oplot, x, hist/double(total(hist)), psym=10,color=l*255d0/ss.nchains;, transparency=100d0-100d0/ss.nchains + + endfor + + hist = histogram(pars,nbins=100,locations=x,min=xmin,max=xmax) + oplot, x, hist/double(total(hist)), psym=10, thick=3 ;; if the best parameters are given, overplot them on the PDFs if finite(best) then $ diff --git a/exofastv2.pro b/exofastv2.pro index 3794762d4..1db6a0225 100644 --- a/exofastv2.pro +++ b/exofastv2.pro @@ -377,7 +377,7 @@ pro exofastv2, priorfile=priorfile, $ maxgr=maxgr, mintz=mintz, $ noyy=noyy, noclaret=noclaret, tides=tides, nplanets=nplanets, $ fitrv=fitrv, fittran=fittran,ttvs=ttvs, earth=earth,$ - i180=i180, covar=covar,alloworbitcrossing=alloworbitcrossing + i180=i180, covar=covar,alloworbitcrossing=alloworbitcrossing,stretch=stretch ;; this is the stellar system structure COMMON chi2_block, ss @@ -410,8 +410,7 @@ for i=0, n_tags(ss)-1 do begin endfor endfor -nchains = n_elements((*ss.tofit)[0,*])*2 -memrequired = nchains*maxsteps*npars*8d0/(1024d0^3) +memrequired = ss.nchains*maxsteps*npars*8d0/(1024d0^3) print, 'Fit will require ' + strtrim(memrequired,2) + ' GB of RAM for the final structure' if memrequired gt 2d0 then begin print, 'WARNING: this likely exceeds your available RAM and may crash after the end of a very long run. You likely want to reduce MAXSTEPS and increase NTHIN by the same factor. If you would like to proceed anyway, type ".con" to continue' @@ -515,9 +514,9 @@ for i=0, ss.nplanets-1 do begin endif if ~ss.planet[i].fitrv and ss.planet[i].fittran then begin ss.planet[i].mpearth.value = massradius_chenreverse(ss.planet[i].rpearth.value) - ss.planet[i].mp.value = ss.planet[i].mpearth.value*0.00000300245 ;; m_sun - ss.planet[i].k.value = (2d0*!dpi*G/(ss.planet[i].period.value*(ss.star.mstar.value + ss.planet[i].mp.value)^2))^(1d0/3d0)*$ - ss.planet[i].mp.value*sin(ss.planet[i].i.value)/sqrt(1d0-ss.planet[i].e.value)*rsun/86400d0 ;; m/s + ss.planet[i].mpsun.value = ss.planet[i].mpearth.value*0.00000300245 ;; m_sun + ss.planet[i].k.value = (2d0*!dpi*G/(ss.planet[i].period.value*(ss.star.mstar.value + ss.planet[i].mpsun.value)^2))^(1d0/3d0)*$ + ss.planet[i].mpsun.value*sin(ss.planet[i].i.value)/sqrt(1d0-ss.planet[i].e.value)*rsun/86400d0 ;; m/s ss.planet[i].logk.value = alog10(ss.planet[i].k.value) endif endif @@ -542,7 +541,7 @@ if not keyword_set(bestonly) then begin nthin=nthin,maxsteps=maxsteps,$ burnndx=burnndx, seed=seed, randomfunc=randomfunc, $ gelmanrubin=gelmanrubin, tz=tz, maxgr=maxgr, mintz=mintz, $ - derived=rstar + derived=rstar,stretch=stretch if pars[0] eq -1 then begin printf, lun, 'MCMC Failed to find a stepping scale. This usually means one or more parameters are unconstrained by the data or priors.' if lun ne -1 then free_lun, lun @@ -589,6 +588,7 @@ mcmcss = mkss(rvpath=rvpath, tranpath=tranpath, fluxfile=fluxfile, nplanets=npla rossiter=rossiter,longcadence=longcadence, earth=earth, i180=i180,$ chen=chen,nvalues=nsteps*nchains,/silent,noyy=noyy,noclaret=noclaret,$ alloworbitcrossing=alloworbitcrossing) +mcmcss.nchains = nchains mcmcss.burnndx = burnndx if ~ss.noyy then mcmcss.star.rstar.value = rstar *(mcmcss.chi2) = chi2 @@ -607,6 +607,7 @@ label = "tab:" + basename caption = "Median values and 68\% confidence interval for " + basename parfile = prefix + 'pdf.ps' covarfile = prefix + 'covar.ps' +chainfile = prefix + 'chain.ps' texfile = prefix + 'median.tex' if keyword_set(covar) then nocovar = 0 $ @@ -614,6 +615,7 @@ else nocovar = 1 exofast_plotdist2, mcmcss, pdfname=parfile, covarname=covarfile,nocovar=nocovar exofast_latextab2, mcmcss, caption=caption, label=label,texfile=texfile +exofast_plotchains, mcmcss, chainfile=chainfile ;; display all the plots, if desired if keyword_set(display) then begin diff --git a/massradius_yy3.pro b/massradius_yy3.pro index 2defe6732..18b79c695 100644 --- a/massradius_yy3.pro +++ b/massradius_yy3.pro @@ -159,11 +159,12 @@ yyrstar = interpol(yytrack[1,*], yytrack[2,*],age, quad=quad, lsquad=lsquad) if yyrstar le 0 then return, !values.d_infinity ;chi2 = ((yyteff-teff)/uteff)^2 -chi2 = ((yyteff-teff)/(yyteff*0.01))^2 +chi2 = ((yyteff-teff)/(yyteff*0.01))^2 ;; this assumes a 1% error in the YY tracks -;; prepare the plotting device -mydevice=!d.name if keyword_set(debug) or keyword_set(psname) then begin + mydevice=!d.name + + ;; prepare the plotting device if keyword_set(psname) then begin ;; astrobetter.com tip on making pretty IDL plots set_plot, 'PS' @@ -179,7 +180,6 @@ if keyword_set(debug) or keyword_set(psname) then begin xtitle=textoidl('T_{eff}') ytitle=textoidl('log g_*') endif else begin -; set_plot, 'X' red = '0000ff'x symsize = 1 device,window_state=win_state @@ -188,35 +188,25 @@ if keyword_set(debug) or keyword_set(psname) then begin xtitle='T_eff' ytitle='log g' endelse -endif - -if keyword_set(debug) or keyword_set(psname) then begin + ;; make a publication-ready plot of the YY track -- Teff vs logg G = 2942.71377d0 ;; R_sun^3/(m_sun*day^2), Torres 2010 rstar = yytrack[1,*] loggplottrack = alog10(G*mstar/(rstar^2)*9.31686171d0) teffplottrack = yytrack[0,*] loggplot = alog10(G*mstar/(yyrstar^2)*9.31686171d0) - xmin=max(teffplottrack,min=xmax) ;; plot range backwards ymin = min([loggplot,3,5],max=ymax) plot, teffplottrack, loggplottrack,xtitle=xtitle,ytitle=ytitle, xrange=[xmin,xmax], yrange=[ymin,ymax] plotsym,0,/fill - oplot, [teff], [loggplot], psym=8,symsize=0.5 - -; plot, alog10(yytrack[0,*]),yytrack[2,*] -; plot, alog10(yyteffall), yyloggall, xrange=[xmax,xmin],yrange=yrange,xtitle=xtitle,ytitle=ytitle -; oploterror, alog10([teff]), [logg], 0.00222,0.01,/lobar -; oploterror, alog10([teff]), [logg], 0.00316978,0.011,/hibar + oplot, [teff], [loggplot], psym=8,symsize=0.5 ;; the model point -endif + if keyword_set(psname) then device, /close + set_plot, mydevice -if keyword_set(psname) then begin - device, /close endif -set_plot, mydevice return, chi2 diff --git a/mkss.pro b/mkss.pro index 270f71eae..6466791cb 100644 --- a/mkss.pro +++ b/mkss.pro @@ -264,6 +264,14 @@ Av.fit = 0 Av.derive = 0 Av.scale = 0.3d0 +alpha = parameter +alpha.description = 'Alpha abundance' +alpha.latex = '\alpha' +alpha.label = 'alpha' +alpha.fit = 0 +alpha.derive = 0 +alpha.scale = 0.3d0 + ;Ma = parameter ;Ma.unit = '' ;Ma.description = 'Apparent V-band Magnitude' @@ -970,6 +978,7 @@ star = create_struct(mstar.label,mstar,$ vsini.label,vsini,$ macturb.label,macturb,$ Av.label,Av,$ + alpha.label,alpha,$ ; Ma.label,Ma,$ ; Mv.label,Mv,$ errscale.label,errscale,$ @@ -1134,6 +1143,7 @@ ss = create_struct('star',star,$ 'ttvs', ttvs,$ 'alloworbitcrossing', alloworbitcrossing,$ 'nsteps',nsteps,$ + 'nchains',1L,$ ; 'nbad',0L,$ 'burnndx',0L,$ 'amoeba',0L,$ @@ -1322,9 +1332,11 @@ if ~keyword_set(silent) then print, 'Those with "no prior constraint" only affec if ~keyword_set(silent) then print openr, lun, priorfile, /get_lun line = '' -i=0 +lineno=0 while not eof(lun) do begin + lineno+=1 readf, lun, line + ;; skip commented lines if strpos(line,'#') eq 0 then continue @@ -1334,7 +1346,7 @@ while not eof(lun) do begin nentries = n_elements(entries) ;; each line must have at least a name and value if nentries lt 2 or nentries gt 5 then begin - message, 'WARNING: line ' + strtrim(i,2) + ' in ' + priorfile + ' is not legal syntax (NAME VALUE [UNCERTAINTY] [LOWERBOUND] [UPPERBOUND]); ignoring: ' + line, /continue + if line ne '' then message, 'WARNING: line ' + strtrim(lineno,2) + ' in ' + priorfile + ' is not legal syntax (NAME VALUE [UNCERTAINTY] [LOWERBOUND] [UPPERBOUND]); ignoring: ' + line, /continue continue endif @@ -1380,15 +1392,21 @@ while not eof(lun) do begin ;; priorwidth = 0 => fix it at the prior value ss.(i)[priornum].(ndx).fit = 0d0 if ~keyword_set(silent) then print, priorname + ' = ' + strtrim(priorval,2) + ' (fixed)' - endif else if finite(priorwidth) and priorwidth gt 0d0 then begin + endif else if finite(priorwidth) and priorwidth gt 0d0 or finite(lowerbound) or finite(upperbound) then begin ;; apply a Gaussian prior with width = priorwidth - ss.(i)[priornum].(ndx).priorwidth = priorwidth priors=[[priors],[i,priornum,ndx]] - ss.(i)[priornum].(ndx).scale = priorwidth*3d0 - + + if priorwidth gt 0 and finite(priorwidth) then begin + ss.(i)[priornum].(ndx).scale = priorwidth*3d0 + ss.(i)[priornum].(ndx).priorwidth = priorwidth + endif + if ~keyword_set(silent) then print, priorname + ' = ' + strtrim(priorval,2) + ' +/- ' + strtrim(priorwidth,2) + '; bounded between ' + strtrim(lowerbound,2) + ' and ' + strtrim(upperbound,2) endif else begin + + + ;; else no prior, just change the default starting value if ~keyword_set(silent) then print, priorname + ' = ' + strtrim(priorval,2) + ' (no prior constraint); bounded between ' + strtrim(lowerbound,2) + ' and ' + strtrim(upperbound,2) endelse @@ -1451,6 +1469,7 @@ for i=0L, n_tags(ss)-1 do begin endfor tofit = tofit[*,1:*] *(ss.tofit) = tofit +ss.nchains = n_elements((*ss.tofit)[0,*])*2L ;; determine the epoch for each observation for i=0, ntran-1 do begin diff --git a/pars2derived.pro b/pars2derived.pro index b96de7f68..4bfa5c1b5 100644 --- a/pars2derived.pro +++ b/pars2derived.pro @@ -18,7 +18,7 @@ for i=0, ss.nplanets-1 do begin ss.planet[i].i.value, ss.planet[i].period.value, $ ss.star.mstar.value) - ss.planet[i].arsun.value=(G*(ss.star.mstar.value+ss.planet[i].mp.value)*ss.planet[i].period.value^2/$ + ss.planet[i].arsun.value=(G*(ss.star.mstar.value+ss.planet[i].mpsun.value)*ss.planet[i].period.value^2/$ (4d0*!dpi^2))^(1d0/3d0) ;; rsun ss.planet[i].ar.value = ss.planet[i].arsun.value/ss.star.rstar.value ;; unitless ss.planet[i].a.value = ss.planet[i].arsun.value/AU ;; AU diff --git a/plotrv.pro b/plotrv.pro index f12ff60d1..82afe51c8 100644 --- a/plotrv.pro +++ b/plotrv.pro @@ -95,9 +95,9 @@ for i=0, ss.nplanets-1 do begin endfor ; xtitle=TeXtoIDL('Phase + (T_P - T_C)/P + 0.25') - xtitle1='Phase + (T_P - T_C)/P + 0.25' + xtitle1='!3Phase + (T_P - T_C)/P + 0.25' plot, [0], [0], xrange=[0,1], yrange=[allminrv,allmaxrv], $ - xtitle=xtitle1,ytitle='RV (m/s)', charsize=charsize, title=ss.planet[i].label + xtitle=xtitle1,ytitle='!3RV (m/s)', charsize=charsize, title=ss.planet[i].label for j=0, ss.ntel-1 do begin rv = *(ss.telescope[j].rvptrs) modelrv = exofast_rv(rv.bjd,ss.planet[i].tp.value,$ @@ -130,10 +130,10 @@ if 0 then begin roundto = 10L^(strlen(strtrim(floor(allmaxdate-allmindate),2))+1L) bjd0 = floor(allmindate/roundto)*roundto -xtitle2='BJD_TDB - ' + string(bjd0,format='(i7)') +xtitle2='!3BJD_TDB - ' + string(bjd0,format='(i7)') plot, [0], [0], xrange=[allmindate,allmaxdate]-bjd0, $ yrange=[allminrv,allmaxrv], xtitle=xtitle2,$ - ytitle='RV (m/s)' + ytitle='!3RV (m/s)' ;TeXtoIDL('BJD_{TDB} - ') + string(bjd0,format='(i7)'), oplot, prettytime-bjd0, allprettymodel, color=red diff --git a/step2pars.pro b/step2pars.pro index 7e1bcddf2..2e937167e 100644 --- a/step2pars.pro +++ b/step2pars.pro @@ -61,8 +61,20 @@ for i=0, ss.nplanets-1 do begin ss.planet[i].e.value = 0d0 ss.planet[i].qesinw.value = 0d0 ss.planet[i].qecosw.value = 0d0 + endif ;else if ss.planet[i].e.value gt (1d0-1d0/ss.planet[i].ar.value) + + ;; limit eccentricity to avoid collision with star during periastron + ;; the ignored tidal effects would become important long before this, + ;; but this prevents numerical problems compared to the e < 1 constraint + ;; the not/lt (instead of ge) robustly handles NaNs too + ;; abs(p) because p is allowed to be negative to eliminate bias + if not (ss.planet[i].e.value lt (1d0-1d0/ss.planet[i].ar.value-abs(ss.planet[i].p.value)/$ + ss.planet[i].ar.value)) then begin + if keyword_set(verbose) then print, 'Planet ' + strtrim(i,2) + ' will collide with the star! e=' + strtrim(ss.planet[i].e.value,2) + '; a/Rstar=' + strtrim(ss.planet[i].ar.value,2) + '; Rp/Rstar=' + strtrim(ss.planet[i].p.value,2) + '; Rstar=' + strtrim(ss.star.rstar.value,2) + if keyword_set(verbose) then print, 'Rstar is derived from YY isochrones; adjust starting values for Age, Teff, [Fe/H], or Mstar to change it' + return, -1 endif - + if ss.planet[i].e.value eq 0d0 then ss.planet[i].omega.value = !dpi/2d0 if ~finite(ss.planet[i].ar.value ) then stop @@ -80,18 +92,6 @@ for i=0, ss.nplanets-1 do begin return, -1 endif - ;; limit eccentricity to avoid collision with star during periastron - ;; the ignored tidal effects would become important long before this, - ;; but this prevents numerical problems compared to the e < 1 constraint - ;; the not/lt (instead of ge) robustly handles NaNs too - ;; abs(p) because p is allowed to be negative to eliminate bias - if not (ss.planet[i].e.value lt (1d0-1d0/ss.planet[i].ar.value-abs(ss.planet[i].p.value)/$ - ss.planet[i].ar.value)) then begin - if keyword_set(verbose) then print, 'Planet ' + strtrim(i,2) + ' will collide with the star! e=' + strtrim(ss.planet[i].e.value,2) + '; a/Rstar=' + strtrim(ss.planet[i].ar.value,2) + '; Rp/Rstar=' + strtrim(ss.planet[i].p.value,2) + '; Rstar=' + strtrim(ss.star.rstar.value,2) - if keyword_set(verbose) then print, 'Rstar is derived from YY isochrones; adjust starting values for Age, Teff, [Fe/H], or Mstar to change it' - return, -1 - endif - ;; for multi-planet systems, make sure they don't enter each other's hill spheres ;; if mp unknown, mp=0 => hill radius=0 => planets can't cross orbits ;; **** ignores mutual inclination; a priori excludes systems like Neptune and Pluto!! **** diff --git a/target2bjd.pro b/target2bjd.pro index 68aeec2bd..e2913b762 100755 --- a/target2bjd.pro +++ b/target2bjd.pro @@ -87,8 +87,8 @@ trueanom = 2d0*atan(sqrt((1d0 + e)/(1d0 - e))*tan(0.5d0*eccanom)) ;; displacement from barycenter if n_elements(q) ne 0 then begin - if keyword_set(primary) then factor = 1d0/(1d0+q) $ - else factor = (1d0 - 1d0/(1d0+q)) + if keyword_set(primary) then factor = 1d0/(1d0+q) $ ;; a*factor = a1 + else factor = q/(1d0+q) ;; a*factor = a2 endif else factor = 1d0 ;; infinite mass ratio ;; distance from barycenter to target