isochroneshttps://github.com/timothydmorton/isochrones
Stellar model grid interpolation and star-property fitting
Use grids of stellar models to...
...conditioned on observations (broadband photometry, spectroscopy, asteroseismology, parallax, etc.)
Precomputed model grids (e.g. the MIST models) contain tables of observable and theoretical properties of stars as a function of various gridded parameters. These grids can be organized in different ways, e.g.
*The MIST and Dartmouth grids define and use a quantity called EEP (equivalent evolutionary phase), which allows for regular gridding in the third dimension, and thus can serve as a proxy for mass (in isochrone grids) or age (in the evolutionary track grids).
Observed properties, e.g.:
$$ \mathbf x = \{\bar J, \bar H, \bar K, \bar \pi\}$$Parameters of model (single star):
$$\mathbf \theta = \{EEP, T, [Fe/H], d, AV \}~~~\mathrm{[using~isochrone~grids]} $$or $$\mathbf \theta = \{M, EEP, [Fe/H], d, AV \}~~~\mathrm{[using~evolution~track~grids]} $$
Incorporate uncertainties in likelihood:
$$p(\mathbf x~|~\mathbf \theta) \propto \prod ...$$Computation of likelihood requires prediction of $x_i$ at arbitrary $\theta$ → Requires interpolation
In [1]:
from isochrones import get_ichrone
mist = get_ichrone('mist', bands='JHK') # Downloads & reorganizes appropriate data into ~/.isochrones (or $ISOCHRONES)
mist.df.head()
Out[1]:
In [2]:
mist.interp
Out[2]:
In [3]:
mist.interp.grid.shape # filled-in regular grid (nan-padded)
Out[3]:
In [4]:
mist.interp.index_columns # feh, log(age), eep
Out[4]:
In [5]:
pars = [0.01, 9.54, 300.3]
mist.interp(pars, 'log_g')
Out[5]:
In [6]:
from scipy.interpolate import RegularGridInterpolator
# Construct SciPy regular grid interpolator for logg
points = mist.interp.index_columns
values = mist.interp.grid[:, :, :, 5] # logg is column 5
fn = RegularGridInterpolator(points, values)
In [7]:
pars = [0.01, 9.54, 300.3]
fn(pars)
Out[7]:
In [8]:
%timeit fn(pars)
In [9]:
mist.interp(pars, 'log_g')
Out[9]:
In [10]:
# convenience function for mist.interp(pars, 'log_g'), with different args
%timeit mist.logg(*pars[::-1])
Slightly slower than SciPy for vectorized calculations, however
In [11]:
import numpy as np
N = 10000
%timeit mist.interp([np.ones(N)*0.01, np.ones(N)*9.54, np.ones(N)*300.3], 'log_g')
In [12]:
%timeit fn(np.array([np.ones(N)*0.01, np.ones(N)*9.54, np.ones(N)*300.3]).T)
In [13]:
from isochrones import StarModel
from isochrones import get_ichrone
mist = get_ichrone('mist', bands='JHK')
props = dict(Teff=(5642, 50.0), feh=(-0.27, 0.08), logg=(4.443, 0.028),
J=(10.523, 0.02), H=(10.211, 0.02), K=(10.152, 0.02))
# Single star model
mod = StarModel(mist, **props)
mod.print_ascii()
In [14]:
mod.param_names
Out[14]:
In [15]:
p = [300, 9.5, 0.1, 300, 0.1]
mod.lnpost(p)
Out[15]:
In [16]:
%timeit mod.lnpost(p)
In [17]:
mod.fit(basename='kep22') # Defaults to using MultiNest if available
In [18]:
mod.samples.describe()
Out[18]:
Generate a corner plot with mod.corner_physical():
In [19]:
mod2 = StarModel(mist, N=2, **props)
mod2.print_ascii()
In [20]:
mod2.param_names
Out[20]:
→ Can also do triple systems, resolved binaries, partially resolved binaries, etc.
In [21]:
%%file demo_star/star.ini
Teff = 4135, 98.0
feh = -0.46, 0.16
logg = 4.711, 0.1
[twomass]
J = 13.513, 0.02
H = 12.845, 0.02
K = 12.693, 0.02
[NIRC2]
resolution = 0.1
separation_1 = 0.6
PA_1 = 100
K_1 = 3.66, 0.05
H_1 = 3.77, 0.03
J_1 = 3.74, 0.05
separation_2 = 1.2
PA_2 = 200
K_2 = 5.1, 0.1
H_2 = 5.2, 0.1
J_2 = 5.15, 0.1
In [22]:
mod3 = StarModel.from_ini(mist, 'demo_star')
mod3.print_ascii()
In [23]:
mod3.param_names
Out[23]:
In [24]:
mod4 = StarModel.from_ini(mist, 'demo_star', N=[2,1,1])
mod4.print_ascii()
In [25]:
mod4.param_names
Out[25]:
In [26]:
mod5 = StarModel.from_ini(mist, 'demo_star', N=[2,1,1], index=[0,1,1])
mod5.print_ascii()
In [27]:
mod5.param_names
Out[27]:
In [28]:
from isochrones import get_ichrone
from isochrones.cluster import SimulatedCluster
mist = get_ichrone('mist')
N = 50
cluster_pars = [8.84, -0.2, 500, 0.03, -3, 0.3, 0.3]
stars = SimulatedCluster(N, *cluster_pars, bands='gri', mass_range=(1, 2.5), phot_unc=0.01)
stars.df.head()
Out[28]:
In [29]:
from isochrones.cluster import StarClusterModel
model = StarClusterModel(mist, stars, eep_bounds=(200, 700), minq=0.5)
model.param_names
Out[29]:
In [30]:
isinstance(model, StarModel)
Out[30]:
In [31]:
model.lnpost(cluster_pars)
Out[31]:
In [32]:
%timeit model.lnpost(cluster_pars)
→ Yikes! What's going on here?
In [33]:
import holoviews as hv
hv.extension('bokeh')
In [34]:
%%opts Points [width=400, height=400, tools=['hover']]
from isochrones.cluster import StarCatalog
import pandas as pd
df = pd.read_hdf('overview/small-test-cluster.h5')
test_stars = StarCatalog(df, bands='gri', props=['parallax'])
test_stars.hr
Out[34]:
In [35]:
%%opts Points [width=500, height=500, tools=['hover']]
data = hv.Points(test_stars.ds, kdims=['g-i', 'g_mag'],
vdims=['is_binary', 'distance',
'mass_pri', 'mass_sec',
'eep_pri', 'eep_sec'], label='Simulated data').options(size=5)
model = mist.hr_isochrone('g', 'i', *cluster_pars[:4], mineep=300, maxeep=700, thin=2, label='Model isochrone').options(size=2)
data * model
Out[35]:
isochrones
In [ ]: