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')

See: Sample K cuts.ipynb for setting StarResult.high_K


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))

Bulk properties


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()

Example light curves


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

High-$K$ examples


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)

Low-$K$ examples


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 [ ]: