In [3]:
from exosyspop.tests.test_populations import test_binpop, test_bgpop, test_planets
from exosyspop.tests.test_populations import get_binpop, get_bgpop, get_planets
In [13]:
import logging
rootLogger = logging.getLogger()
rootLogger.setLevel(logging.DEBUG)
In [51]:
#pop = get_bgpop()
pop = get_planets()
In [67]:
from exosyspop.populations import PoissonPlanetPopulation
import cPickle as pickle
rB = pop.radius_B
s1 = pop.stars
rootLogger.setLevel(logging.INFO)
pop._train_trap(N=200)
rootLogger.setLevel(logging.DEBUG)
obs = pop.observe(regr_trap=True)
pop.save('test', overwrite=True)
#pickle.dump(pop, open('test.pkl', 'wb'))
#pop = pickle.load(open('test.pkl', 'rb'))
pop = PoissonPlanetPopulation('test')
s2 = pop.stars
rB2 = pop.radius_B
In [73]:
(rB==rB2).all()
Out[73]:
In [69]:
cols = ['mass_A', 'mass_B', 'radius_A', 'radius_B','period']
s1[cols].head(15)
Out[69]:
In [70]:
s2[cols].head(15)
Out[70]:
In [56]:
(s1 == s2).all()
Out[56]:
In [20]:
len(rB2)
Out[20]:
In [22]:
sum(rB==rB2)/float(len(rB))
Out[22]:
In [ ]:
In [16]:
rB[:10]
Out[16]:
In [17]:
rB2[:10]
Out[17]:
In [6]:
test_binpop()
test_bgpop()
test_planets()
In [1]:
import sys
## This to get the peerless target star DataFrame for example purposes
#sys.path.append('/u/tdm/repositories/peerless/prediction')
#sys.path.append('/u/tdm/repositories/peerless')
#from targets import targets
import pandas as pd
targets = pd.read_hdf('targets.h5')
# Sanitize dtypes of targets DataFrame
for c in targets.columns:
if targets[c].dtype == object:
targets.loc[:,c] = targets.loc[:,c].astype(str)
# The action is here. Depends on vespa & isochrones.
from exosyspop.populations import KeplerBinaryPopulation
In [2]:
#from isochrones.dartmouth import Dartmouth_Isochrone
pop = KeplerBinaryPopulation(targets, fB=0.4)
In [3]:
pop._train_trap(N=200)
Out[3]:
In [3]:
from exosyspop.tests.test_populations import _do_the_things
In [4]:
_do_the_things(pop)
In [9]:
# Accessing secondary properties will initialize a secondary simulation,
# calling pop._generate_binaries(). The first time this is called, the
# secondary property regressors get trained.
pop.radius_B
Out[9]:
In [11]:
pop._train_trap(N=200)
Out[11]:
In [12]:
obs = pop.observe(regr_trap=True)
In [13]:
pop.save('testpop', overwrite=True)
In [14]:
import cPickle as pickle
pickle.dump(pop, open('testpop.pkl','wb'))
In [15]:
pop = pickle.load(open('testpop.pkl','rb'))
In [16]:
obs = pop.observe(regr_trap=True)
In [7]:
pop.radius_B
Out[7]:
In [17]:
pop.__dict__.keys()
Out[17]:
In [4]:
# subsequent calls are much faster; e.g.
pop._generate_binaries()
print(pop.radius_B)
%timeit pop._generate_binaries()
In [5]:
# If physical accuracy is important, you can also choose to generate binary properties
# directly from the isochrone, but it's a factor of a few slower:
pop._generate_binaries(use_ic=True)
print(pop.radius_B)
%timeit pop._generate_binaries(use_ic=True)
In [7]:
# Similarly, accessing orbital properties will generate them
pop.period
Out[7]:
In [9]:
# Now, we can observe and see what we see. This takes into account
# duty cycle & data span, as well as geometry.
obs = pop.observe()
print(len(obs))
print(obs.columns)
obs.head()
Out[9]:
In [10]:
# This is pretty fast, even when generating a new population each time:
%timeit pop.observe(new=True)
In [11]:
# Even faster if we only generate new orbits.
%timeit pop.observe(new_orbits=True)
In [12]:
# So we can predict the expected number of observations pretty easily.
import numpy as np
N = 100
n_obs = np.array([len(pop.observe(new_orbits=True)) for i in range(N)])
n_obs.mean(), n_obs.std()
Out[12]:
In [13]:
# Notice that the above does not yet have trapezoidal parameters. There are two options to generate these.
# Either we can set the fit_trap parameter, as follows:
obs = pop.observe(fit_trap=True)
print(len(obs))
obs.columns
Out[13]:
In [14]:
# All things considered, this is still pretty fast if we just need to do it a few times:
%timeit pop.observe(fit_trap=True)
In [15]:
# However, this is pretty slow if we want to do inference. To help with this, we can
# tell it to train & use a regression. Training only happens once; by default with 10,000
# synthetic observations. This takes a minute or so.
obs = pop.observe(regr_trap=True)
print(len(obs))
obs.columns
Out[15]:
In [16]:
# Subsequent calls are much faster
%timeit pop.observe(regr_trap=True)
In [17]:
# Even generating a new stellar population & observing it is pretty quick
%timeit pop.observe(regr_trap=True, new=True)
In [16]:
# Or again, you can just generate new orbits (rather than new binaries & new orbits)
%timeit pop.observe(regr_trap=True, new_orbits=True)
In [40]:
# Generating the training data used for the trapezoid shape regression above used
# this function, which can be otherwise useful to sample >N random observations
# from the existing population. `trap_regr` defaults to `True` here.
# This function also takes `new` or `new_orbits` keywords.
obs_pop = pop.get_N_observed(N=2000, new_orbits=True)
print(len(obs_pop))
obs_pop.columns
Out[40]:
In [45]:
# Let's look at the distribution of trapezoid parameters
%matplotlib inline
from corner import corner
def trap_corner(df, regr=True, **kwargs):
pri = df.n_pri.values > 0
sec = (df.n_sec.values > 0) & ~pri
if regr:
dur = pri*df.trap_dur_pri_regr.fillna(0) + sec*df.trap_dur_sec_regr.fillna(0)
logd = np.log10(pri*df.trap_depth_pri_regr.fillna(0) +
sec*df.trap_depth_sec_regr.fillna(0))
slope = pri*df.trap_slope_pri_regr.fillna(0) + sec*df.trap_slope_sec_regr.fillna(0)
else:
dur = pri*df.trap_dur_pri.fillna(0) + sec*df.trap_dur_sec.fillna(0)
logd = np.log10(pri*df.trap_depth_pri.fillna(0) + sec*df.trap_depth_sec.fillna(0))
slope = pri*df.trap_slope_pri.fillna(0) + sec*df.trap_slope_sec.fillna(0)
#return (dur, logd, slope)
return corner(np.array([dur, logd, slope]).T, labels=['duration', 'log(d)', 'T/tau'],
plot_contours=False, **kwargs)
trap_corner(obs_pop, range=[(0,0.6),(-4,0), (2,6)]);
In [42]:
obs_pop2 = pop.get_N_observed(N=2000, new_orbits=True, fit_trap=True)
In [46]:
trap_corner(obs_pop2, regr=False, range=[(0,0.6),(-4,0), (2,6)]);
In [18]:
# We can now look, e.g. at the expected number of single/double eclipsing systems:
query = '(n_pri < 3) & (n_sec < 3) & (n_pri==0 | n_sec==0)'
N = 100
n_obs = np.array([len(pop.observe(new_orbits=True).query(query)) for i in range(N)])
n_obs.mean(), n_obs.std()
Out[18]:
In [19]:
# Try this again, this time using the empirical eccentricity distribution
# (as opposed to the beta distribution with default params)---eccentricity matters!
pop.ecc_empirical = True
n_obs = np.array([len(pop.observe(new_orbits=True).query(query)) for i in range(N)])
n_obs.mean(), n_obs.std()
Out[19]:
In [20]:
# You can also save a trained model so that you can load it back and hit the ground running.
pop.save('ebpop', overwrite=True)
pop = KeplerBinaryPopulation.load('ebpop')
In [21]:
# No training necessary!
pop.observe(regr_trap=True).head()
Out[21]:
In [ ]: