<jupyter_start><jupyter_text>trace plots<jupyter_code>ax = []
fig = plt.figure(figsize=(20,12))

for pp, par in enumerate(par_bwm):
    ax.append(fig.add_subplot(2,3,pp+1))
    ax[pp].plot(chain_bwm[:,pp])
    ax[pp].set_xlabel(par)

ax.append(fig.add_subplot(2,3,6))
ax[5].plot(chain[:,-4])
ax[5].set_xlabel('log_likelihood');<jupyter_output><empty_output><jupyter_text>Corner and Histograms<jupyter_code>corner_kwargs = {'bins':30,
                 'show_titles':True,
                 'labels':par_bwm,
                 'smooth':1.0,
                 'plot_datapoints':False,
                 'plot_density':True,
                 'plot_contours':True,
                 'fill_contours':False,}

fig = corner(chain_bwm, color='C0', **corner_kwargs);
ax = []
fig = plt.figure(figsize=(20,8))

for pp, par in enumerate(par_bwm):
    ax.append(fig.add_subplot(2,3,pp+1))
    ax[pp].hist(chain_bwm[:,pp], bins=25, density=True, histtype='step') #, label='')
    ax[pp].set_xlabel(par)

#ax[0].legend(loc='upper left')
ax[idx_A].set_yscale('log')
ax[idx_t0].set_yscale('log')

ax.append(fig.add_subplot(2,3,6))
ax[5].hist(chain[:,-4], bins=30, density=True, histtype='step')
ax[5].set_xlabel('log_likelihood')
lognorm = LogNorm()

hist_kwargs = {'bins':30,
               'normed':True,
               'cmap':'Blues',
               #'vmin':1.0e5,
               #'vmax':1.0e10,
               'norm':lognorm,
                }

fig = plt.figure(figsize=(8,4))
ax = fig.add_subplot(111)

t_ch = chain_bwm[:,idx_t0]
logA_ch = chain_bwm[:,idx_A]
A_ch = 10**logA_ch
    
ax.hist2d(t_ch, logA_ch, **hist_kwargs)

if ii_t:
    ax.axvline(x=tmin, color='k', linestyle=':')
    ax.axvline(x=tmax, color='k', linestyle=':')

ax.set_xlabel(r'$t_0$ (MJD)')
ax.set_ylabel(r'$\log_{10}(A_{BWM})$')
ax.set_ylim([-18,-11])
ax.set_xlim([TMIN, TMAX])
#ax.set_title('{:.1f} yr timeslice (w/ dropout)'.format(slice_yr))
import matplotlib.gridspec as gridspec
from matplotlib.ticker import NullFormatter

UL, dUL = UL_uncert(A_ch, p=0.95)

# Define size of figure
fig = plt.figure(figsize=(20, 15))
gs = gridspec.GridSpec(10, 12)

# Define the positions of the subplots.
ax0 = plt.subplot(gs[6:10, 5:9])
axx = plt.subplot(gs[5:6, 5:9])
axy = plt.subplot(gs[6:10, 9:10])


hist_kwargs = {'bins':30,
               'normed':True,
               'cmap':'Blues',
              }
ax0.hist2d(t_ch, A_ch, **hist_kwargs)
ax0.set_xlabel(r'$t_0$ (MJD)')
ax0.set_ylabel(r'$A_{BWM}$')
ax0.set_xlim([TMIN, TMAX])
#ax0.set_ylim([0, 0.13e-12])

axx.hist(t_ch, 30, normed=True, histtype='step', color='C0')
axx.set_yscale('log')
axy.hist(A_ch, 30, normed=True, histtype='step', color='C0', orientation='horizontal')
axy.axhline(y=UL, color='C1', linestyle='--')
axy.set_xscale('log')
               
# Remove tick labels
nullfmt = NullFormatter()
axx.xaxis.set_major_formatter(nullfmt)
axx.yaxis.set_major_formatter(nullfmt)
axy.xaxis.set_major_formatter(nullfmt)
axy.yaxis.set_major_formatter(nullfmt)<jupyter_output>/home/pbaker/.local/opt/miniconda2/lib/python2.7/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.
  warnings.warn("The 'normed' kwarg is deprecated, and has been "<jupyter_text>UL v. t_0<jupyter_code>burnfrac = 0.15
thin = 20

ephems = ['DE421', 'DE430', 'DE436']
bayeseph_flags = [True, False]

#ephems = ['DE421']
#bayeseph_flags = [False]
#figs = []

ULs = {}
dULs = {}

for ephem in ephems:
    for bayesephem in bayeseph_flags:
        BE = '_BE' if bayesephem else ''
        key = ephem+BE
        print(key)
        ULs[key] = []
        dULs[key] = []

        for ts in tlim:
            tmin, cent, tmax = ts
            chunk = '{:.2f}'.format(cent)
            rundir = runroot + '{0:s}/{1:s}/'.format(key, chunk)
            chain_raw = pd.read_csv(rundir + 'chain_1.txt',
                                    sep='\t', dtype=float, header=None).values
            burn = int(burnfrac * len(chain_raw))
            chain = chain_raw[burn::thin]
            chain_bwm = chain[:,idx_bwm]  # no ephem model
            chain_L = chain[:,-4]

            corL = acor(chain_L)[0]
            corA = acor(chain_bwm[:,idx_A])[0]
            N = len(chain_bwm)
            print("{:.2f}: N = {:d}, corL = {:.2f}, corA = {:.2f}".format(cent, N, corL, corA))
    
            this_UL, this_dUL = UL_uncert(10**chain_bwm[:,idx_A], p=0.95)
            ULs[key].append(this_UL)
            dULs[key].append(this_dUL)
            
            #figs.append(hist_plot(chain_bwm, par_bwm, cols=5, wid_per_col=2.5))

print('DONE!')
# save per chunk hists
outdir = "hists_"+key
!mkdir -p $outdir

for ii,fig in enumerate(figs):
    fig.savefig("hists_"+key+"/hist_{:02d}.png".format(ii))
    plt.close(fig)
plt.rcParams.update(plt.rcParamsDefault)
params = {'backend': 'pdf',
        'axes.labelsize': 10,
        'lines.markersize': 4,
        'font.size': 10,
        'xtick.major.size':6,
        'xtick.minor.size':3,
        'ytick.major.size':6,
        'ytick.minor.size':3,
        'xtick.major.width':0.7,
        'ytick.major.width':0.7,
        'xtick.minor.width':0.7,
        'ytick.minor.width':0.7,
        'lines.markeredgewidth':1.2,
        'axes.linewidth':1.375,
        'legend.fontsize': 7,
        'xtick.labelsize': 12,
        'ytick.labelsize': 12,
        'savefig.dpi':200,
        'path.simplify':True,
        #'font.family': 'serif',
        #'font.serif':'Times',
        'text.latex.preamble':'\\usepackage{amsmath}',
        'text.usetex':True,
         }

plt.rcParams.update(params)
from astropy.time import Time

Ts = np.array(tlim)[:,1]

Ts = Time(Ts, format='mjd').decimalyear
plt.close('all')

fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111)

ax.errorbar(Ts, ULs['DE421'], yerr=dULs['DE421'], color='C0', label='DE421')
ax.errorbar(Ts, ULs['DE430'], yerr=dULs['DE430'], color='C1', label='DE430')
ax.errorbar(Ts, ULs['DE436'], yerr=dULs['DE436'], color='C2', label='DE436')

#ax.errorbar(Ts, ULs['DE421_BE'], yerr=dULs['DE421_BE'], color='C0', label='DE421')
#ax.errorbar(Ts, ULs['DE430_BE'], yerr=dULs['DE430_BE'], color='C1', label='DE430')
#ax.errorbar(Ts, ULs['DE436_BE'], yerr=dULs['DE436_BE'], color='C1', label='DE436')

ax.legend(loc='lower left', fontsize=16)
ax.set_yscale('log')
ax.set_xlabel(r'burst epoch', fontsize=16)
ax.set_ylabel(r'$95\%$ UL (strain)', fontsize=16)

fig.savefig('/home/pbaker/11yr_ULvt.pdf')
plt.show()
plt.close('all')

fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111)

ax.errorbar(Ts, ULs['DE421'], yerr=dULs['DE421'], color='C0', alpha=0.5, ls='--')
ax.errorbar(Ts, ULs['DE430'], yerr=dULs['DE430'], color='C1', alpha=0.5, ls='--')
ax.errorbar(Ts, ULs['DE436'], yerr=dULs['DE436'], color='C2', alpha=0.5, ls='--')

ax.errorbar(Ts, ULs['DE421_BE'], yerr=dULs['DE421_BE'], color='C0', label='DE421')
ax.errorbar(Ts, ULs['DE430_BE'], yerr=dULs['DE430_BE'], color='C1', label='DE430')
ax.errorbar(Ts, ULs['DE436_BE'], yerr=dULs['DE436_BE'], color='C2', label='DE436')

ax.legend(loc='lower left', fontsize=16)
ax.set_yscale('log')
ax.set_xlabel(r'burst epoch', fontsize=16)
ax.set_ylabel(r'$95\%$ UL (strain)', fontsize=16)

fig.savefig('/home/pbaker/11yr_ULvt_BE.pdf')
plt.show()
plt.close('all')

ephem = 'DE436'

fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111)

ax.errorbar(Ts, ULs[ephem], yerr=dULs[ephem], color='C0', alpha=0.5, ls='--')
ax.errorbar(Ts, ULs[ephem+'_BE'], yerr=dULs[ephem+'_BE'], color='C0', label=ephem)

ax.legend(loc='lower left', fontsize=16)
ax.set_yscale('log')
ax.set_xlabel(r'burst epoch', fontsize=16)
ax.set_ylabel(r'$95\%$ UL (strain)', fontsize=16)

#fig.savefig('/home/pbaker/11yr_ULvt_compare.pdf')
plt.show()
plt.close('all')

fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111)

ax.plot(Ts, ULs['DE421'], color='C0', alpha=0.5, ls='--')
ax.plot(Ts, ULs['DE430'], color='C1', alpha=0.5, ls='--')
ax.plot(Ts, ULs['DE436'], color='C2', alpha=0.5, ls='--')

ax.plot(Ts, ULs['DE421_BE'], color='C0', label='DE421')
ax.plot(Ts, ULs['DE430_BE'], color='C1', label='DE430')
ax.plot(Ts, ULs['DE436_BE'], color='C2', label='DE436')

ax.legend(loc='lower left', fontsize=16)
ax.set_yscale('log')
ax.set_xlabel(r'burst epoch', fontsize=16)
ax.set_ylabel(r'$95\%$ UL (strain)', fontsize=16)

#fig.savefig('/home/pbaker/11yr_ULvt_BE.pdf')
plt.show()<jupyter_output><empty_output><jupyter_text>Rate Upper Limit<jupyter_code>from scipy.integrate import simps

_GMsun = 1.32712440018e20  # m^3/s^2
_c = 299792458 # m/s
_Rsun = _GMsun / _c**2 # meters
_pc = 3.0857e16 # meters

M = 1e10 * _Rsun / (1e9*_pc) # Gpc
Erad = 0.07*M  # energy radiated

# h = 69/(8*24) Erad/r

hUL = np.array(ULs['DE436_BE'])

r_sens = 69/(8*24) * Erad / hUL
V_sens = 4/3*np.pi * r_sens**3

rateUL = simps(y=V_sens, x=Ts) # /Gpc^3 /yr
69/(8*24) * 0.07 *_Rsun
rateUL
plt.plot(Ts, r_sens)
dT = np.mean(np.diff(Ts))
hsort = np.sort(hUL)
Ns = np.arange(len(hsort))+1
plt.plot(hsort, 1/(Ns*dT))
plt.xlabel('strain, $h_\mathrm{gw}$')
plt.ylabel('event rate UL (N/yr)')
plt.xscale('log')
plt.yscale('log')
cols = [k for k in ULs.keys()]
ULvt = np.vstack([u for u in ULs.values()])
dULvt = np.vstack([u for u in dULs.values()])

np.savetxt('./ULvT/9yr_cols.dat', cols, fmt=b'%s')
np.savetxt('./ULvT/9yr_UL.dat', ULvt)
np.savetxt('./ULvT/9yr_dUL.dat', dULvt)
np.savetxt('./ULvT/9yr_Ts.dat', Ts)
Tobs = TMAX-TMIN

0.1*Tobs / 365.25
UL<jupyter_output><empty_output>