In [1]:
import holoviews as hv
hv.extension('bokeh')
In [2]:
# Simulate a cluster
import pickle
import numpy as np
import pandas as pd
from isochrones import get_ichrone
from isochrones.priors import PowerLawPrior
from isochrones.utils import addmags
from isochrones.cluster import simulate_cluster
resim = True
if resim:
N = 50
age = 8.0
feh = -0.42
distance = 200.
AV = 0.0
alpha = -2
gamma = 0.3
fB = 0.5
cat = simulate_cluster(N, age, feh, distance, AV, alpha, gamma, fB)
cat.df.to_hdf('test_cluster.h5', 'df')
pardict = dict(age=age, feh=feh, distance=distance, AV=AV, alpha=alpha, gamma=gamma, fB=fB)
pickle.dump(pardict, open('test_cluster_params.pkl', 'wb'))
else:
bands = 'JHK'
stars = pd.read_hdf('test_cluster.h5')
pardict = pickle.load(open('test_cluster_params.pkl', 'rb'))
age = pardict['age']
feh = pardict['feh']
distance = pardict['distance']
AV = pardict['AV']
alpha = pardict['alpha']
gamma = pardict['gamma']
fB = pardict['fB']
In [3]:
cat.df.head()
Out[3]:
In [4]:
import holoviews as hv
hv.extension('bokeh')
cat.hr
Out[4]:
In [5]:
import pandas as pd
from isochrones.cluster import StarClusterModel
from isochrones import get_ichrone
from isochrones.priors import FehPrior
from isochrones.mist import MIST_Isochrone
mist = MIST_Isochrone()
model = StarClusterModel(mist, cat, eep_bounds=(202, 605),
max_distance=3000, max_AV=0.1, name='test-binary')
model.set_prior('feh', FehPrior(halo_fraction=0.5))
print(model.param_names)
In [6]:
print(model.mnest_basename)
In [7]:
pars = [age, feh, distance, AV, alpha, gamma, fB]
model.lnprior(pars), model.lnlike(pars), model.lnpost(pars)
In [10]:
model.ic.interp_mag([[300, 301, 302, 303], 9.6, 0.0, 200, 0.1], bands=['J', 'H', 'K'])
Out[10]:
In [11]:
cat.df.head()
Out[11]:
In [ ]:
In [13]:
from isochrones import StarModel
from isochrones.mist import MIST_Isochrone, MIST_EvolutionTrack
iso = MIST_Isochrone(bands=['PS_g', 'PS_r', 'PS_i', 'PS_z', 'PS_y', 'PS_w'])
mod = StarModel(iso, J=(8.43, 0.02), H=(8.32, 0.02), K=(8.30, 0.02), parallax=(4.92, 0.1))
In [43]:
iso.initialize() # Downloads/unpacks data, builds local grids, etc.
In [14]:
mod.param_names
Out[14]:
In [15]:
mod2 = StarModel(iso, J=(8.43, 0.02), H=(8.32, 0.02), K=(8.30, 0.02), parallax=(4.92, 0.1), N=2)
mod2.param_names
Out[15]:
In [22]:
pars = [300, 9.6, 0.0, 200, 0.1]
mod.lnlike(pars), mod.lnprior(pars), mod.lnpost(pars)
Out[22]:
In [25]:
pars2 = [330, 300, 9.7, 0.0, 200, 0.2]
mod2.lnpost(pars2)
Out[25]:
In [27]:
mod.fit(basename='mytest')
In [32]:
%matplotlib inline
import matplotlib as mpl
mpl.style.use('dark_background')
mod.corner_physical();
In [33]:
mod.prior_transform?
In [35]:
cube = [0.2] * 5
mod.prior_transform(cube)
Out[35]:
In [38]:
from isochrones.mist import MIST_EvolutionTrack
from isochrones.starmodel import BasicStarModel
track = MIST_EvolutionTrack()
mod_evtrack = BasicStarModel(track, J=(8.43, 0.02), H=(8.32, 0.02), K=(8.30, 0.02), parallax=(4.92, 0.1))
In [39]:
mod_evtrack.param_names
Out[39]:
In [41]:
pars_evtrack = [1.2, 300, 0.0, 200, 0.1]
mod_evtrack.lnlike(pars_evtrack)
Out[41]:
In [42]:
%timeit mod_evtrack.lnlike(pars_evtrack)
In [ ]:
pars2 = []
In [9]:
pars2 = [9.5, -0.52, 200., 0.02, -1.5, 0.3, 0.6]
model.lnprior(pars2), model.lnlike(pars2), model.lnpost(pars2)
Out[9]:
In [8]:
cube = [0.5] * 7
model.mnest_prior(cube, 7, 7)
model.mnest_loglike(cube, 7, 7)
Out[8]:
In [9]:
pickle.dump(model, open('testmodel.pkl', 'wb'))
model = pickle.load(open('testmodel.pkl', 'rb'))
In [10]:
model.lnpost(pars)
Out[10]:
In [26]:
model.ic._interp
Out[26]:
In [8]:
%timeit model.lnpost([age, feh, distance, AV, alpha, gamma, 0.5])
In [15]:
model._priors
Out[15]:
In [7]:
nsamples = 20
for i in range(nsamples):
# pardict = {k : pr.sample(1)[0] for k,pr in model._priors.items()}
# pars = [pardict[k] for k in model.param_names]
cube = [np.random.random() for p in model.param_names]
model.mnest_prior(cube, 7, 7)
print('{}: {}'.format([float('{:.2f}'.format(x)) for x in cube], model.mnest_loglike(cube, 7, 7)))
In [9]:
rand_pars = [9.5, 0.0, 500, 0.1, -2, 0.5, 0.4]
model.lnpost(rand_pars)
Out[9]:
In [9]:
import pickle
pickle.dump(model, open('testmodel.pkl', 'wb'))
In [10]:
model = pickle.load(open('testmodel.pkl', 'rb'))
In [11]:
model.lnpost(pars)
In [ ]:
model.fit(overwrite=True)
In [ ]:
! say "sampling complete!"
In [30]:
%matplotlib inline
import corner
names = ['age', 'feh', 'distance', 'AV', 'gamma']
fig = corner.corner(model.samples[names], names=names,
truths=[age, feh, distance, AV, gamma])
In [ ]:
In [13]:
bands = 'rJK'
def band_pairs(bands):
return [(bands[i], bands[-1]) for i in range(len(bands)-1)]
band_pairs(bands)
Out[13]:
In [14]:
def iso_compare(ic, bands, age1, feh1, dist1, AV1,
age2, feh2, dist2, AV2):
def make_df(iso):
df = pd.DataFrame()
for b1, b2 in band_pairs(bands):
mag1 = iso['{}_mag'.format(b1)]
mag2 = iso['{}_mag'.format(b2)]
df[b2] = mag2
df['{0}-{1}'.format(b1, b2)] = mag1 - mag2
return df
iso1 = ic.isochrone(age1, feh1, distance=dist1, AV=AV1)
ds1 = hv.Dataset(make_df(iso1))
iso2 = ic.isochrone(age2, feh2, distance=dist2, AV=AV2)
ds2 = hv.Dataset(make_df(iso2))
layout = []
for b1, b2 in band_pairs(bands):
kdims = ['{}-{}'.format(b1, b2), b2]
layout.append(hv.Points(ds1, kdims) * hv.Points(ds2, kdims))
return hv.Layout(layout)
In [15]:
%%opts Points {+framewise}
from functools import partial
dmap = hv.DynamicMap(partial(iso_compare, ic=mist, bands='rJK'),
kdims=['age1', 'feh1', 'dist1', 'AV1', 'age2', 'feh2', 'dist2', 'AV2'])
age_bounds = model.bounds('age')
feh_bounds = model.bounds('feh')
dist_bounds = model.bounds('distance')
dmap.redim.range(age1=age_bounds, age2=age_bounds,
feh1=feh_bounds, feh2=feh_bounds,
dist1=dist_bounds, dist2=dist_bounds)
Out[15]:
In [16]:
means = model.samples.mean()
iso_compare(mist, bands='rJK', age1=age, feh1=feh, dist1=distance, AV1=AV,
age2=means.age, feh2=means.feh, dist2=means.distance, AV2=means.AV)
Out[16]:
In [18]:
age2, feh2, dist2, AV2, gamma2 = other_pars
In [19]:
iso_compare(mist, bands='rJK', age1=age, feh1=feh, dist1=distance, AV1=AV,
age2=age2, feh2=feh2, dist2=dist2, AV2=AV2)
Out[19]:
In [15]:
means = model.samples.mean()
means
Out[15]:
In [16]:
mean_pars = [6.85, -1.0, 13300, 0.01, -0.94]
mean_pars = [means.age, means.feh, means.distance, means.AV, means.gamma]
model.lnprior(mean_pars), model.lnlike(mean_pars), model.lnpost(mean_pars)
Out[16]:
In [17]:
model.lnprior(pars), model.lnlike(pars), model.lnpost(pars)
Out[17]:
In [19]:
%debug model.lnlike(pars)
In [ ]: