<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Occurrence Rate Calculations (Binned)
Step1:  Planet Database
Step2:  G Dwarfs
Load up the G-star table
Step3:  We want to know the average detection efficiency across the bin.  (Actually, what we want to know is the fraction of the planets appearing in the bin that would be detected, but since we are assuming that the planet distribution function is flat, this reduces to the average.)
$$
\alpha_i \equiv \frac{1}{\Delta P \Delta R} \int dP\,dR\, p_{det}(P,R)
$$
For each of the Kep ID's in the G-star table, load up the efficiency contour, and calculate the efficiency integrated over the relevant bin limits.  Save a two-column kepid and alpha file.  There are lots of errors because there are a number of stars that pass our cuts that do not have efficiency calculations from Chris---I think these are the un-classified stars that have $T$ and $log(g)$ replaced with solar values.
Step4:  Now load the data file we just generated
Step5:  The designated bin is the (2,2) bin (1.5 to 2.25 in $R$, 40-80 in $P$).  Here is a histogram of the average detection efficiency in the designated bin
Step6:  The posterior on the rate is proportional to
$$
p(R | N_{pl}, \alpha) \sim R^{N_{pl}} \exp\left[ - R \sum_i \alpha_i \right]
$$
Let's grab the number of counts in the bins
Step7:  M-Dwarfs
Repeat the above for the M-dwarfs
Step8:  We are interested in the (1, 2) bin for M-dwarfs.  What does the detection efficiency look like there?
Step9:  And Finally for K-Dwarfs
<ASSISTANT_TASK:>
Python Code:
hw_data_directory = '/Users/farr/Documents/Research/KepHackWeek/data'
occur_dir = '/Users/farr/Google Drive/Kepler ExoPop Hack 2015/end2end_occ_calc'
eff_dir = '/Volumes/KepHacWkWMF/Kepler_HW2015/Dp4_DetectionCountours/v0'
rbins = array([1.5**(i-1) for i in range(9)])
pbins = array([10*2**i for i in range(6)])
print rbins, pbins
def integrate_efficiency(ps, rs, eff, pmin, pmax, rmin, rmax):
    sel = (ps >= pmin) & (ps < pmax)
    p_integral = sint.simps(eff[sel, :], ps[sel], axis=0)
    p_integral = p_integral / (np.max(ps[sel])-np.min(ps[sel]))
    sel = (rs >= rmin) & (rs < rmax)
    return sint.simps(p_integral[sel], rs[sel]) / (np.max(rs[sel]) - np.min(rs[sel]))
def rate_and_uncertainty(npl, nexpected):
    Rbest = npl / nexpected
    dr = sqrt(npl + 1) / nexpected
    
    Rmin = max(0.0, Rbest - 10*dr)
    Rmax = Rbest + 10*dr
    
    Rs = linspace(Rmin, Rmax, 10000)
    logpost = (npl-0.5)*log(Rs) - nexpected*Rs
    post = exp(logpost)
    post = post / trapz(post, Rs)
    cumpost = sint.cumtrapz(post, Rs)
    
    imin = argmin(abs(cumpost - 0.16))
    imed = argmin(abs(cumpost - 0.5))
    imax = argmin(abs(cumpost - 0.84))
    
    return Rs[imed], Rs[imin], Rs[imax]
Explanation: Occurrence Rate Calculations (Binned)
End of explanation
planets = aa.read(op.join(hw_data_directory, 'dr24_koi.csv'))
planets
Explanation: Planet Database
End of explanation
gdata = genfromtxt(op.join(occur_dir, 'hack_week_g_stars.csv'), delimiter=',', names=True)
Explanation: G Dwarfs
Load up the G-star table:
End of explanation
with gzip.open('g-integrated-eff.dat.gz.temp', 'w') as out:
    out.write('# kepid ')
    out.write(' '.join(['(%g-%g, %g-%g)'%(plow, phigh, rlow, rhigh) \
                        for (plow, phigh) in zip(pbins[:-1], pbins[1:]) \
                        for (rlow, rhigh) in zip(rbins[:-1], rbins[1:])]))
    out.write('\n')
    for kid in gdata['kepid']:
        try:
            kid = int(round(kid))
            eff = af.open(op.join(eff_dir, 'contourV0%09d_tot.fits.gz'%kid))
            ps = exp(linspace(log(eff[0].header['MINPER']), log(eff[0].header['MAXPER']), eff[0].header['NPER']))
            rs = exp(linspace(log(eff[0].header['MINRP']), log(eff[0].header['MAXRP']), eff[0].header['NRP']))
            int_effs = []
            for plow, phigh in zip(pbins[:-1], pbins[1:]):
                for rlow, rhigh in zip(rbins[:-1], rbins[1:]):
                    int_effs.append(integrate_efficiency(ps, rs, eff[0].data.T, plow, phigh, rlow, rhigh))
            int_effs = array(int_effs)
            out.write('%09d '%kid)
            out.write(' '.join(['%g'%ie for ie in int_effs]))
            out.write('\n')
            eff.close()
        except:
            print 'Could not process KID ', kid, ' exception: ', sys.exc_info()[0]
os.rename('g-integrated-eff.dat.gz.temp', 'g-integrated-eff.dat.gz')
Explanation: We want to know the average detection efficiency across the bin.  (Actually, what we want to know is the fraction of the planets appearing in the bin that would be detected, but since we are assuming that the planet distribution function is flat, this reduces to the average.)
$$
\alpha_i \equiv \frac{1}{\Delta P \Delta R} \int dP\,dR\, p_{det}(P,R)
$$
For each of the Kep ID's in the G-star table, load up the efficiency contour, and calculate the efficiency integrated over the relevant bin limits.  Save a two-column kepid and alpha file.  There are lots of errors because there are a number of stars that pass our cuts that do not have efficiency calculations from Chris---I think these are the un-classified stars that have $T$ and $log(g)$ replaced with solar values.
End of explanation
g_det_eff = loadtxt('g-integrated-eff.dat.gz')
g_det_ids = g_det_eff[:,0]
g_det_eff = g_det_eff[:,1:].reshape((-1, pbins.shape[0]-1, rbins.shape[0]-1))
Explanation: Now load the data file we just generated:
End of explanation
pu.plot_histogram_posterior(g_det_eff[:,2,2])
xlabel(r'$\alpha$')
ylabel(r'$p(\alpha)$')
Explanation: The designated bin is the (2,2) bin (1.5 to 2.25 in $R$, 40-80 in $P$).  Here is a histogram of the average detection efficiency in the designated bin:
End of explanation
gkids_set = set(g_det_ids)
g_pl_counts = zeros((pbins.shape[0]-1, rbins.shape[0]-1), dtype=np.int)
noutside = 0
for pl in planets:
    if pl['kepid'] in gkids_set and pl['koi_pdisposition'] == 'CANDIDATE' and pl['koi_max_mult_ev'] >= 15:
        pind = bisect.bisect(pbins, pl['koi_period']) - 1
        rind = bisect.bisect(rbins, pl['koi_prad']) - 1
                
        if (pind < 0) or (rind < 0) or (pind >= g_pl_counts.shape[0]) or (rind >= g_pl_counts.shape[1]):
            noutside += 1
        else:
            g_pl_counts[pind, rind] += 1
print noutside, ' planets landed outside any bins!'
ggammas = zeros((pbins.shape[0]-1, rbins.shape[0]-1))
ggammamins = zeros((pbins.shape[0]-1, rbins.shape[0]-1))
ggammamaxs = zeros((pbins.shape[0]-1, rbins.shape[0]-1))
for i in range(pbins.shape[0]-1):
    for j in range(rbins.shape[0]-1):
        rmed, rmin, rmax = rate_and_uncertainty(g_pl_counts[i,j], np.sum(g_det_eff[:,i,j]))
        ggammas[i,j] = rmed
        ggammamins[i,j] = rmin
        ggammamaxs[i,j] = rmax
PS, RS = meshgrid(pbins, rbins, indexing='ij')
pcolormesh(PS, RS, ggammas, norm=LogNorm())
xscale('log')
yscale('log')
colorbar()
axis(xmin=min(pbins), xmax=max(pbins), ymin=min(rbins), ymax=max(rbins))
xlabel(r'$P$ (d)')
ylabel(r'$R$ ($R_\oplus$)')
savefig('eta_g.pdf')
with gzip.open('pbins.dat.gz', 'w') as out:
    savetxt(out, pbins.reshape((1,-1)))
with gzip.open('rbins.dat.gz', 'w') as out:
    savetxt(out, rbins.reshape((1, -1)))
with gzip.open('eta_g.dat.gz', 'w') as out:
    savetxt(out, ggammas)
with gzip.open('sigma_p_g.dat.gz', 'w') as out:
    savetxt(out, ggammamaxs-ggammas)
with gzip.open('sigma_n_g.dat.gz', 'w') as out:
    savetxt(out, ggammas-ggammamins)
Explanation: The posterior on the rate is proportional to
$$
p(R | N_{pl}, \alpha) \sim R^{N_{pl}} \exp\left[ - R \sum_i \alpha_i \right]
$$
Let's grab the number of counts in the bins:
End of explanation
mdata = genfromtxt(op.join(occur_dir, 'hack_week_m_dwarfs.csv'), delimiter=',', names=True)
with gzip.open('m-integrated-eff.dat.gz.temp', 'w') as out:
    out.write('# kepid ')
    out.write(' '.join(['(%g-%g, %g-%g)'%(plow, phigh, rlow, rhigh) \
                        for (plow, phigh) in zip(pbins[:-1], pbins[1:]) \
                        for (rlow, rhigh) in zip(rbins[:-1], rbins[1:])]))
    out.write('\n')
    for kid in mdata['kepid']:
        try:
            kid = int(round(kid))
            eff = af.open(op.join(eff_dir, 'contourV0%09d_tot.fits.gz'%kid))
            ps = exp(linspace(log(eff[0].header['MINPER']), log(eff[0].header['MAXPER']), eff[0].header['NPER']))
            rs = exp(linspace(log(eff[0].header['MINRP']), log(eff[0].header['MAXRP']), eff[0].header['NRP']))
            int_effs = []
            for plow, phigh in zip(pbins[:-1], pbins[1:]):
                for rlow, rhigh in zip(rbins[:-1], rbins[1:]):
                    int_effs.append(integrate_efficiency(ps, rs, eff[0].data.T, plow, phigh, rlow, rhigh))
            int_effs = array(int_effs)
            out.write('%09d '%kid)
            out.write(' '.join(['%g'%ie for ie in int_effs]))
            out.write('\n')
            eff.close()
        except:
            print 'Could not process KID ', kid, ' exception: ', sys.exc_info()[0]
os.rename('m-integrated-eff.dat.gz.temp', 'm-integrated-eff.dat.gz')
m_det_eff = loadtxt('m-integrated-eff.dat.gz')
m_det_ids = m_det_eff[:,0]
m_det_eff = m_det_eff[:,1:].reshape((-1, pbins.shape[0]-1, rbins.shape[0]-1))
Explanation: M-Dwarfs
Repeat the above for the M-dwarfs
End of explanation
pu.plot_histogram_posterior(m_det_eff[:,1,2])
mkids_set = set(m_det_ids)
m_pl_counts = zeros((pbins.shape[0]-1, rbins.shape[0]-1), dtype=np.int)
noutside = 0
for pl in planets:
    if pl['kepid'] in mkids_set and pl['koi_pdisposition'] == 'CANDIDATE' and pl['koi_max_mult_ev'] >= 15:
        pind = bisect.bisect(pbins, pl['koi_period']) - 1
        rind = bisect.bisect(rbins, pl['koi_prad']) - 1
                
        if (pind < 0) or (rind < 0) or (pind >= m_pl_counts.shape[0]) or (rind >= m_pl_counts.shape[1]):
            noutside += 1
        else:
            m_pl_counts[pind, rind] += 1
print noutside, ' planets landed outside any bins!'
mgammas = zeros((pbins.shape[0]-1, rbins.shape[0]-1))
mgammamins = zeros((pbins.shape[0]-1, rbins.shape[0]-1))
mgammamaxs = zeros((pbins.shape[0]-1, rbins.shape[0]-1))
for i in range(pbins.shape[0]-1):
    for j in range(rbins.shape[0]-1):
        rmed, rmin, rmax = rate_and_uncertainty(m_pl_counts[i,j], np.sum(m_det_eff[:,i,j]))
        mgammas[i,j] = rmed
        mgammamins[i,j] = rmin
        mgammamaxs[i,j] = rmax
PS, RS = meshgrid(pbins, rbins, indexing='ij')
pcolormesh(PS, RS, mgammas, norm=LogNorm())
xscale('log')
yscale('log')
colorbar()
axis(xmin=min(pbins), xmax=max(pbins), ymin=min(rbins), ymax=max(rbins))
xlabel(r'$P$ (d)')
ylabel(r'$R$ ($R_\oplus$)')
savefig('eta_m.pdf')
print 'In the (1,2) bin R = ', mgammas[1,2], ' + ', (mgammamaxs[1,2]-mgammas[1,2]), ' - ', (mgammas[1,2] - mgammamins[1,2])
with gzip.open('pbins.dat.gz', 'w') as out:
    savetxt(out, pbins.reshape((1,-1)))
with gzip.open('rbins.dat.gz', 'w') as out:
    savetxt(out, rbins.reshape((1, -1)))
with gzip.open('eta_m.dat.gz', 'w') as out:
    savetxt(out, mgammas)
with gzip.open('sigma_p_m.dat.gz', 'w') as out:
    savetxt(out, mgammamaxs-mgammas)
with gzip.open('sigma_n_m.dat.gz', 'w') as out:
    savetxt(out, mgammas-mgammamins)
Explanation: We are interested in the (1, 2) bin for M-dwarfs.  What does the detection efficiency look like there?
End of explanation
kdata = genfromtxt(op.join(occur_dir, 'hack_week_k_dwarfs.csv'), delimiter=',', names=True)
with gzip.open('k-integrated-eff.dat.gz.temp', 'w') as out:
    out.write('# kepid ')
    out.write(' '.join(['(%g-%g, %g-%g)'%(plow, phigh, rlow, rhigh) \
                        for (plow, phigh) in zip(pbins[:-1], pbins[1:]) \
                        for (rlow, rhigh) in zip(rbins[:-1], rbins[1:])]))
    out.write('\n')
    for kid in kdata['kepid']:
        try:
            kid = int(round(kid))
            eff = af.open(op.join(eff_dir, 'contourV0%09d_tot.fits.gz'%kid))
            ps = exp(linspace(log(eff[0].header['MINPER']), log(eff[0].header['MAXPER']), eff[0].header['NPER']))
            rs = exp(linspace(log(eff[0].header['MINRP']), log(eff[0].header['MAXRP']), eff[0].header['NRP']))
            int_effs = []
            for plow, phigh in zip(pbins[:-1], pbins[1:]):
                for rlow, rhigh in zip(rbins[:-1], rbins[1:]):
                    int_effs.append(integrate_efficiency(ps, rs, eff[0].data.T, plow, phigh, rlow, rhigh))
            int_effs = array(int_effs)
            out.write('%09d '%kid)
            out.write(' '.join(['%g'%ie for ie in int_effs]))
            out.write('\n')
            eff.close()
        except:
            print 'Could not process KID ', kid, ' exception: ', sys.exc_info()[0]
os.rename('k-integrated-eff.dat.gz.temp', 'k-integrated-eff.dat.gz')
k_det_eff = loadtxt('k-integrated-eff.dat.gz')
k_det_ids = k_det_eff[:,0]
k_det_eff = k_det_eff[:,1:].reshape((-1, pbins.shape[0]-1, rbins.shape[0]-1))
kkids_set = set(k_det_ids)
k_pl_counts = zeros((pbins.shape[0]-1, rbins.shape[0]-1), dtype=np.int)
noutside = 0
for pl in planets:
    if pl['kepid'] in kkids_set and pl['koi_pdisposition'] == 'CANDIDATE' and pl['koi_max_mult_ev'] >= 15:
        pind = bisect.bisect(pbins, pl['koi_period']) - 1
        rind = bisect.bisect(rbins, pl['koi_prad']) - 1
                
        if (pind < 0) or (rind < 0) or (pind >= k_pl_counts.shape[0]) or (rind >= k_pl_counts.shape[1]):
            noutside += 1
        else:
            k_pl_counts[pind, rind] += 1
print noutside, ' planets landed outside any bins!'
kgammas = zeros((pbins.shape[0]-1, rbins.shape[0]-1))
kgammamins = zeros((pbins.shape[0]-1, rbins.shape[0]-1))
kgammamaxs = zeros((pbins.shape[0]-1, rbins.shape[0]-1))
for i in range(pbins.shape[0]-1):
    for j in range(rbins.shape[0]-1):
        rmed, rmin, rmax = rate_and_uncertainty(k_pl_counts[i,j], np.sum(k_det_eff[:,i,j]))
        kgammas[i,j] = rmed
        kgammamins[i,j] = rmin
        kgammamaxs[i,j] = rmax
PS, RS = meshgrid(pbins, rbins, indexing='ij')
pcolormesh(PS, RS, kgammas, norm=LogNorm())
xscale('log')
yscale('log')
colorbar()
axis(xmin=min(pbins), xmax=max(pbins), ymin=min(rbins), ymax=max(rbins))
xlabel(r'$P$ (d)')
ylabel(r'$R$ ($R_\oplus$)')
savefig('eta_k.pdf')
print 'In K-star bin (40-80, 1.5-2.25)', kgammas[2, 2], ' + ', (kgammamaxs[2,2]-kgammas[2,2]), ' - ', (kgammas[2,2]-kgammamins[2,2])
print 'In K-star bin (20-40, 1.5-2.25)', kgammas[1, 2], ' + ', (kgammamaxs[1,2]-kgammas[1,2]), ' - ', (kgammas[1,2]-kgammamins[1,2])
with gzip.open('pbins.dat.gz', 'w') as out:
    savetxt(out, pbins.reshape((1,-1)))
with gzip.open('rbins.dat.gz', 'w') as out:
    savetxt(out, rbins.reshape((1, -1)))
with gzip.open('eta_k.dat.gz', 'w') as out:
    savetxt(out, kgammas)
with gzip.open('sigma_p_k.dat.gz', 'w') as out:
    savetxt(out, kgammamaxs-kgammas)
with gzip.open('sigma_n_k.dat.gz', 'w') as out:
    savetxt(out, kgammas-kgammamins)
Explanation: And Finally for K-Dwarfs
End of explanation
<END_TASK>