In [ ]:
import os
from os import path
# Third-party
import astropy.units as u
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import numpy as np
%matplotlib inline
import h5py
import pandas as pd
from sqlalchemy import func
import tqdm
from thejoker import JokerSamples
from twoface.config import TWOFACE_CACHE_PATH
from twoface.samples_analysis import unimodal_P
from twoface.db import (db_connect, AllStar, AllVisit, AllVisitToAllStar,
StarResult, Status, JokerRun, NessRG)
from twoface.plot import plot_two_panel
In [ ]:
plot_path = '../../paper/1-catalog/figures/'
table_path = '../../paper/1-catalog/tables/'
os.makedirs(plot_path, exist_ok=True)
os.makedirs(table_path, exist_ok=True)
Session, _ = db_connect(path.join(TWOFACE_CACHE_PATH, 'apogee.sqlite'))
session = Session()
samples_file = path.join(TWOFACE_CACHE_PATH, 'apogee-jitter.hdf5')
In [ ]:
high_K_stars = session.query(AllStar).join(StarResult).filter(StarResult.status_id>0).filter(StarResult.high_K).all()
low_K_stars = session.query(AllStar).join(StarResult).filter(StarResult.status_id>0).filter(~StarResult.high_K).all()
apogee_ids = np.array([x.apogee_id for x in high_K_stars])
len(high_K_stars), len(low_K_stars)
In [ ]:
n_martig = session.query(AllStar).join(StarResult).filter(StarResult.high_K).filter(AllStar.martig_filter).count()
n_ness = session.query(AllStar).join(StarResult, NessRG).filter(StarResult.high_K).count()
print(n_martig, n_ness)
In [ ]:
statuses = np.array([star.results[0].status_id for star in high_K_stars])
for i in range(4+1):
count = session.query(StarResult).filter(StarResult.status_id == i).count()
status = session.query(Status).filter(Status.id == i).one()
print("Status: {1} [{0}]\n\t{2} in high K sample of {3} total\n"
.format(i, status.message, np.sum(statuses == i), count))
In [ ]:
sp = np.array(session.query(AllStar.logg, AllStar.logg_err, AllStar.teff, AllStar.teff_err)
.join(StarResult).filter(StarResult.status_id>0)
.filter(StarResult.high_K).all())
In [ ]:
fig, ax = plt.subplots(1, 1, figsize=(6, 6))
ax.errorbar(sp[:, 2], sp[:, 0],
xerr=sp[:, 3], yerr=sp[:, 1],
alpha=0.2, marker='o', linestyle='none')
ax.set_xlim(5750, 3450)
ax.set_ylim(4, 0)
In [ ]:
Ks = np.array(session.query(AllStar.k)
.join(StarResult).filter(StarResult.status_id>0)
.filter(StarResult.high_K).all())[:,0]
plt.hist(Ks, bins='auto');
In [ ]:
nvisits = np.array(session.query(func.count(AllVisit.id)).join(AllVisitToAllStar, AllStar, StarResult)
.filter(StarResult.status_id>0)
.filter(StarResult.high_K)
.group_by(AllStar.apogee_id)
.having(func.count(AllVisit.id) >= 3).all())
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
_ = ax.hist(nvisits, bins=np.logspace(0.5, 7, 22, base=2.), rasterized=True)
ax.set_xscale('log', basex=2)
ax.set_yscale('log')
ax.set_xlim(2.5, 140)
ax.set_ylim(0.8, 10**5.1)
ax.xaxis.set_ticks(2**np.arange(2, 7+1, 1))
ax.xaxis.set_ticklabels([str(x) for x in ax.get_xticks()])
ax.yaxis.set_ticks(10**np.arange(0, 5+1, 1))
ax.set_xlabel('$N$ visits')
ax.set_ylabel('$N$ stars')
fig.tight_layout()
In [ ]:
def make_plots(q, title, seed=420):
np.random.seed(seed)
ranges = [np.arange(3, 9+1, 2),
np.arange(11, 14+1, 1)]
rc = {
'axes.labelsize': 18,
'xtick.labelsize': 14,
'ytick.labelsize': 14
}
figs = []
with mpl.rc_context(rc):
for k, range_ in enumerate(ranges):
gs = GridSpec(4, 3)
fig = plt.figure(figsize=(8., 9.5))
for j, i in enumerate(range_):
ax1 = fig.add_subplot(gs[j, :2])
ax2 = fig.add_subplot(gs[j, 2])
if j == 0:
ax1.set_title(title, fontsize=20)
sub_stars = q.having(func.count(AllVisit.id) == i).all()
print(len(sub_stars))
star = np.random.choice(sub_stars)
data = star.apogeervdata()
with h5py.File(samples_file, 'r') as f:
samples = JokerSamples.from_hdf5(f[star.apogee_id])
fig = plot_two_panel(data, samples, axes=[ax1, ax2], tight=False,
plot_data_orbits_kw=dict(n_times=16384,
n_orbits=128,
highlight_P_extrema=False,
xlim_choice='data',
relative_to_t0=True,
plot_kwargs=dict(alpha=0.2,
linewidth=0.2,
rasterized=True)))
xlim = ax1.get_xlim()
ylim = ax1.get_ylim()
ax1.text(xlim[0] + (xlim[1]-xlim[0])/20,
ylim[1] - (ylim[1]-ylim[0])/20,
star.apogee_id, fontsize=15, va='top', ha='left')
ax1.text(xlim[1] - (xlim[1]-xlim[0])/20,
ylim[1] - (ylim[1]-ylim[0])/20,
'$N = {0}$'.format(len(data)),
fontsize=15, va='top', ha='right')
ax1.set_xlabel('')
ax2.set_xlabel('')
ax1.set_xlabel(r'${\rm BMJD} - t_0$ [day]')
ax2.set_xlabel('period, $P$ [day]')
fig.tight_layout()
fig.subplots_adjust(left=0.125, right=0.95, hspace=0.2, wspace=0.4)
figs.append(fig)
return figs
In [ ]:
q = session.query(AllStar).join(StarResult, AllVisitToAllStar, AllVisit)\
.filter(StarResult.status_id>0)\
.filter(StarResult.high_K)\
.filter(StarResult.status_id == 4)\
.filter(AllStar.logg > 1.5)\
.filter(AllStar.aspcapflag.op('&')(2**7) == 0)\
.group_by(AllStar.apstar_id)
figs = make_plots(q, title='Examples of high-$K$ stars')
for i, fig in enumerate(figs):
fig.savefig(path.join(plot_path, 'highK-{0}.pdf'.format(i)),
dpi=150, rasterized=True)
In [ ]:
q = session.query(AllStar).join(StarResult, AllVisitToAllStar, AllVisit)\
.filter(StarResult.status_id>0)\
.filter(~StarResult.high_K)\
.filter(StarResult.status_id == 4)\
.filter(AllStar.logg > 1.5)\
.filter(AllStar.aspcapflag.op('&')(2**7) == 0)\
.group_by(AllStar.apstar_id)
figs = make_plots(q, title='Examples of low-$K$ stars')
for i, fig in enumerate(figs):
fig.savefig(path.join(plot_path, 'lowK-{0}.pdf'.format(i)),
dpi=150, rasterized=True)
In [ ]: