In [1]:
from isochrones.dartmouth import Dartmouth_Isochrone
dar = Dartmouth_Isochrone()
In [2]:
%%file example.obs
name band resolution mag e_mag separation pa relative
twomass J 2.8 10 0.02 0 0 False
twomass H 2.9 9.6 0.02 0 0 False
twomass K 3.0 9.4 0.02 0 0 False
UKIRT J 1.0 11 0.02 0 0 True
UKIRT J 1.0 14.5 0.02 2. 280 True
nirc2 K 0.1 0 0 0 0 True
nirc2 K 0.1 4.0 0.03 2.1 274 True
nirc2 K 0.1 1.5 0.03 0.3 123 True
In [3]:
import pandas as pd
df = pd.read_table('example.obs', delim_whitespace=True)#, index_col=[0,1])
df
Out[3]:
In [4]:
from isochrones.observation import ObservationTree
tree = ObservationTree.from_df(df)
tree.print_ascii()
In [5]:
tree.define_models(dar, N=(2,1,1), index=(0,0,1))
tree.print_ascii()
In [6]:
pars = [1.0,0.9,0.8,9.5,0.0,200,0.1,0.6,9.3,0.1,300,0.2]
pdict = tree.p2pardict(pars)
In [7]:
tree.print_ascii(pars)
In [8]:
tree.add_spectroscopy(Teff=(5700,50), logg=(4.4,0.2))
In [9]:
tree.print_ascii(pars)
In [25]:
tree.lnlike(pars)
Out[25]:
In [26]:
tree.add_spectroscopy(Teff=(5700,50), logg=(4.5,0.2))
In [27]:
tree.lnlike(pars)
Out[27]:
In [32]:
tree.add_parallax((4.8,0.5))
In [33]:
tree.lnlike(pars)
Out[33]:
In [16]:
tree.parallax
Out[16]:
In [15]:
for s,(val,err) in tree.parallax.items():
print s,(val,err)
In [8]:
%timeit tree.get_leaf('0_0')
In [16]:
%timeit tree.get_leaf('0_0').evaluate([1.,9.8,0.0,200,0.3], 'Teff')
In [15]:
%timeit dar.Teff(1,9.8,0.0)
In [13]:
tree.leaf_labels
Out[13]:
In [17]:
tree.systems
Out[17]:
In [7]:
%timeit tree.lnlike(pars)
In [9]:
%timeit dar.mag['K'](1,9.5,0.0,200,0.2)
In [1]:
from isochrones.dartmouth import Dartmouth_Isochrone
from isochrones.starmodel_new import StarModel
import logging
import os
rootLogger = logging.getLogger()
rootLogger.setLevel(logging.INFO)
In [2]:
kep22_props = dict(maxAV = 0.187,
g = (12.0428791, 0.05),
r = (11.5968652, 0.05),
i = (11.4300704, 0.05),
z = (11.393061, 0.05),
J = (10.523, 0.02),
H = (10.211, 0.02),
K = (10.152, 0.02),
Kepler = 11.664,
Teff = (5642, 50.0),
feh = (-0.27, 0.08),
logg = (4.443, 0.028))
mod = StarModel(Dartmouth_Isochrone, N=2, **kep22_props)
In [3]:
mod.obs.print_ascii()
In [4]:
mod.obs.to_df()
Out[4]:
In [7]:
mod.obs.Nstars
Out[7]:
In [5]:
mod.obs.spectroscopy
Out[5]:
In [6]:
mod.obs.parallax
Out[6]:
In [11]:
for o in mod.obs._observations:
for s in o.sources:
print o.name, o.band, o.resolution, s.mag, s.e_mag, s.separation, s.pa, s.relative
In [ ]:
#mod.fit_mcmc()
mod.fit_multinest()
#os.system('say "Fitting is complete"')
In [ ]:
mod.samples
In [4]:
a = mod.mnest_analyzer
In [5]:
a.get_equal_weighted_posterior().shape
Out[5]:
In [6]:
%matplotlib inline
import corner
corner.corner(a.get_equal_weighted_posterior()[:,:-1], labels=['mass_A','mass_B','mass_C','age','feh','dist','AV']);
In [29]:
import matplotlib.pyplot as plt
(a.get_equal_weighted_posterior()[:,2] > a.get_equal_weighted_posterior()[:,1]).sum()
Out[29]:
In [19]:
mod.bounds('mass')
Out[19]:
In [21]:
mod.obs.leaf_labels
Out[21]:
In [6]:
mod.samples.columns
Out[6]:
In [7]:
%matplotlib inline
import matplotlib.pyplot as plt
In [8]:
plt.hist(mod.sampler.flatchain[:,0], bins=50);
In [9]:
plt.hist(mod.sampler.flatchain[:,1], bins=50);
In [11]:
plt.hist(mod.sampler.flatchain[:,-2], bins=50);
In [ ]:
In [22]:
ok = mod.sampler.acceptance_fraction > 0.15
In [24]:
chains = mod.sampler.chain[ok,:,:]
flatchain = chains.reshape((chains.shape[0]*chains.shape[1], chains.shape[2]))
In [30]:
flatchain.shape
Out[30]:
In [15]:
import corner
In [18]:
corner.corner(mod.sampler.flatchain, labels=['mass_A','mass_B','age','feh','dist','AV'],
range=[0.99]*6);
In [27]:
corner?
In [11]:
import numpy as np
w = np.where(mod.sampler.acceptance_fraction < 0.02)[0]
In [14]:
w
Out[14]:
In [13]:
mod.sampler.lnprobability[w,:]
Out[13]:
In [32]:
mod.obs.print_ascii(p=[ 1.92184216e+00, 9.33893417e+00, -1.64027746e+00,
7.95152588e+02, 3.49565065e-02])
In [31]:
mod._sampler.chain[w,:,:]
Out[31]:
In [9]:
mod.obs.print_ascii()
In [10]:
p = [1.0,9.7,0.0,400,0.1]
#p = [1.0,0.5,9.7,0.0,400,0.1]
#p = [1.0,9.7,0.0,400,0.1,0.5,9.5,0.1,200,0.05]
mod.lnpost(p)
Out[10]:
In [12]:
mod.obs.print_ascii(p=p)
In [13]:
mod.
In [5]:
mod.lnprior(p)
Out[5]:
In [7]:
mod.prior('feh',0.0)
Out[7]:
In [8]:
mod.prior('distance',400)
Out[8]:
In [9]:
mod.prior('AV',0.1)
Out[9]:
In [6]:
mod.prior('age',9.7)
Out[6]:
In [6]:
mod.prior('mass',1.0)
Out[6]:
In [12]:
mod._bounds
Out[12]:
In [8]:
for prop in ['age','feh','distance','AV']:
print mod.bounds(prop)
In [20]:
mod.obs.print_ascii(p=p)
In [7]:
import numpy as np
np.log(0)
Out[7]:
In [8]:
np.isinf(-np.inf)
Out[8]:
In [17]:
l = [3,1,20,5,10,2]
m = l[1:4]
m.sort()
m
l
Out[17]:
In [18]:
l
Out[18]:
In [19]:
sorted(l)
Out[19]:
In [11]:
import corner
In [ ]:
corner.