isochrones
¶import holoviews as hv
hv.extension('bokeh')
# 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']
cat.df.head()
J_mag | H_mag | K_mag | is_binary | age | feh | distance | AV | mass_pri | mass_sec | eep_pri | eep_sec | J_mag_unc | H_mag_unc | K_mag_unc | parallax | parallax_unc | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 9.853543 | 9.638057 | 9.604427 | False | 8.0 | -0.42 | 206.100120 | 0.0 | 1.119560 | 0.000000 | 250.709990 | NaN | 0.01 | 0.01 | 0.01 | 4.852011 | 0.2 |
1 | 10.439053 | 10.139190 | 10.117781 | False | 8.0 | -0.42 | 199.473230 | 0.0 | 0.944224 | 0.000000 | 235.834845 | NaN | 0.01 | 0.01 | 0.01 | 5.013204 | 0.2 |
2 | 9.650236 | 9.411557 | 9.405093 | False | 8.0 | -0.42 | 187.979106 | 0.0 | 1.121915 | 0.000000 | 250.925022 | NaN | 0.01 | 0.01 | 0.01 | 5.319740 | 0.2 |
3 | 9.382083 | 9.251201 | 9.252969 | False | 8.0 | -0.42 | 204.870083 | 0.0 | 1.272563 | 0.000000 | 257.903030 | NaN | 0.01 | 0.01 | 0.01 | 4.881142 | 0.2 |
4 | 9.320257 | 9.168760 | 9.143769 | True | 8.0 | -0.42 | 202.396197 | 0.0 | 1.269720 | 0.481506 | 257.687433 | 195.40378 | 0.01 | 0.01 | 0.01 | 4.940804 | 0.2 |
import holoviews as hv
hv.extension('bokeh')
cat.hr
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)
['age', 'feh', 'distance', 'AV', 'alpha', 'gamma', 'fB']
print(model.mnest_basename)
./chains/mist-cluster_test-binary-
pars = [age, feh, distance, AV, alpha, gamma, fB]
model.lnprior(pars), model.lnlike(pars), model.lnpost(pars)
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) <ipython-input-7-aaa46a1915be> in <module>() 1 pars = [age, feh, distance, AV, alpha, gamma, fB] ----> 2 model.lnprior(pars), model.lnlike(pars), model.lnpost(pars) /Users/tdm/repositories/isochrones-dev/isochrones/cluster.py in lnlike(self, p) 373 374 # This takes the lnlike_prop and adds likelihoods from photometry and mass --> 375 mass_lo, mass_hi = self.bounds('mass') 376 Neep = len(eeps) 377 Nbands = len(self.stars.bands) /Users/tdm/repositories/isochrones-dev/isochrones/cluster.py in bounds(self, prop) 276 self.ic.maxeep) 277 elif prop == 'mass': --> 278 return self._mass_bounds if self._mass_bounds is not None else (self.ic.minmass, 279 self.ic.maxmass) 280 AttributeError: 'MIST_Isochrone' object has no attribute 'minmass'
model.ic.interp_mag([[300, 301, 302, 303], 9.6, 0.0, 200, 0.1], bands=['J', 'H', 'K'])
(array([4024.42973009, 4054.35853018, 4085.14927372, 4120.6750109 ]), array([4.70301435, 4.69674803, 4.69214031, 4.68689162]), array([0.02242706, 0.02194302, 0.02160682, 0.02122628]), array([[12.47218125, 11.74176959, 11.58963706], [12.4174364 , 11.68811655, 11.54015964], [12.36842216, 11.64055886, 11.49681227], [12.31250421, 11.58641791, 11.44749982]]))
cat.df.head()
J_mag | H_mag | K_mag | is_binary | age | feh | distance | AV | mass_pri | mass_sec | eep_pri | eep_sec | J_mag_unc | H_mag_unc | K_mag_unc | parallax | parallax_unc | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 8.435677 | 8.325925 | 8.304757 | True | 8.0 | -0.42 | 203.176590 | 0.0 | 1.359026 | 1.356795 | 261.671923 | 261.559643 | 0.01 | 0.01 | 0.01 | 4.921827 | 0.2 |
1 | 9.050041 | 8.961375 | 8.952934 | False | 8.0 | -0.42 | 198.009479 | 0.0 | 1.415914 | 0.000000 | 264.731003 | NaN | 0.01 | 0.01 | 0.01 | 5.050263 | 0.2 |
2 | 8.947427 | 8.893473 | 8.892698 | False | 8.0 | -0.42 | 202.756550 | 0.0 | 1.514407 | 0.000000 | 270.269062 | NaN | 0.01 | 0.01 | 0.01 | 4.932023 | 0.2 |
3 | 9.525714 | 9.348124 | 9.322268 | False | 8.0 | -0.42 | 203.067171 | 0.0 | 1.218387 | 0.000000 | 255.649035 | NaN | 0.01 | 0.01 | 0.01 | 4.924479 | 0.2 |
4 | 8.395321 | 8.408289 | 8.399057 | False | 8.0 | -0.42 | 208.977836 | 0.0 | 2.015060 | 0.000000 | 296.021877 | NaN | 0.01 | 0.01 | 0.01 | 4.785196 | 0.2 |
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))
iso.initialize() # Downloads/unpacks data, builds local grids, etc.
mod.param_names
['eep_0_0', 'age_0', 'feh_0', 'distance_0', 'AV_0']
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
['eep_0_0', 'eep_0_1', 'age_0', 'feh_0', 'distance_0', 'AV_0']
pars = [300, 9.6, 0.0, 200, 0.1]
mod.lnlike(pars), mod.lnprior(pars), mod.lnpost(pars)
(-48587.13049531381, -21.769815906588878, -48608.9003112204)
pars2 = [330, 300, 9.7, 0.0, 200, 0.2]
mod2.lnpost(pars2)
-20186.94986138102
mod.fit(basename='mytest')
INFO:root:MultiNest basename: ./chains/mytest
%matplotlib inline
import matplotlib as mpl
mpl.style.use('dark_background')
mod.corner_physical();
mod.prior_transform?
Signature: mod.prior_transform(cube) Docstring: <no docstring> File: ~/repositories/isochrones-dev/isochrones/starmodel.py Type: method
cube = [0.2] * 5
mod.prior_transform(cube)
array([ 3.420e+02, 6.026e+00, -3.100e+00, 2.000e+03, 2.000e-01])
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))
mod_evtrack.param_names
('mass', 'eep', 'feh', 'distance', 'AV')
pars_evtrack = [1.2, 300, 0.0, 200, 0.1]
mod_evtrack.lnlike(pars_evtrack)
-4860.480635238232
%timeit mod_evtrack.lnlike(pars_evtrack)
1000 loops, best of 3: 233 µs per loop
pars2 = []
pars2 = [9.5, -0.52, 200., 0.02, -1.5, 0.3, 0.6]
model.lnprior(pars2), model.lnlike(pars2), model.lnpost(pars2)
(-12.537630162989952, -601.8697889229954, -614.4074190859853)
cube = [0.5] * 7
model.mnest_prior(cube, 7, 7)
model.mnest_loglike(cube, 7, 7)
-11840.735123021514
pickle.dump(model, open('testmodel.pkl', 'wb'))
model = pickle.load(open('testmodel.pkl', 'rb'))
model.lnpost(pars)
-216.85345003562443
model.ic._interp
<isochrones.interp.DFInterpolator at 0x1c34d00748>
%timeit model.lnpost([age, feh, distance, AV, alpha, gamma, 0.5])
1 loop, best of 3: 385 ms per loop
model._priors
{'age': <isochrones.priors.FlatLogPrior at 0x102852cf8>, 'feh': <isochrones.priors.FehPrior at 0x102852ba8>, 'AV': <isochrones.priors.FlatPrior at 0x1c1bc65080>, 'distance': <isochrones.priors.PowerLawPrior at 0x1c1bc65048>, 'alpha': <isochrones.priors.FlatPrior at 0x1c1bc65320>, 'gamma': <isochrones.priors.FlatPrior at 0x1c1bc652e8>, 'fB': <isochrones.priors.FlatPrior at 0x1c1bc65278>}
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)))
[7.49, -1.36, 840.85, 0.08, -2.61, 0.94, 0.24]: -8650.472829973574 [8.37, -1.89, 1160.27, 0.06, -1.9, 0.25, 0.08]: -10494.460596329816 [6.65, 0.14, 2646.17, 0.04, -2.71, 0.95, 0.48]: -inf [7.0, -2.54, 2028.74, 0.02, -1.06, 0.58, 0.7]: -14140.558981077918 [6.17, -3.35, 2077.46, 0.01, -2.1, 0.99, 0.44]: -14858.678167224747 [7.55, -3.9, 1732.96, 0.1, -3.56, 0.24, 0.42]: -13734.915826425537 [9.0, -2.1, 710.78, 0.04, -3.25, 0.99, 0.88]: -7515.566409828148 [8.57, -0.55, 518.2, 0.07, -2.48, 0.7, 0.22]: -5357.477644060242 [6.99, -3.79, 2740.38, 0.07, -3.0, 0.25, 0.96]: -15948.92165818422 [10.03, -0.82, 2465.53, 0.06, -3.35, 0.26, 0.75]: -10944.405488103577 [7.83, -1.64, 1530.8, 0.07, -2.53, 0.77, 0.98]: -12054.05209673006 [6.42, -2.68, 1694.24, 0.08, -2.39, 0.26, 0.92]: -13369.626149827485 [6.28, -1.17, 2643.34, 0.09, -1.13, 0.65, 0.81]: -14204.942266149212 [10.12, -2.21, 2714.6, 0.07, -3.98, 0.02, 0.23]: -11237.548779934263 [6.35, -2.12, 1713.66, 0.03, -3.12, 0.41, 0.27]: -13459.108961364422 [8.04, -3.77, 960.51, 0.08, -1.93, 0.02, 0.01]: -9924.362713057426 [9.07, -1.41, 2481.95, 0.01, -2.01, 0.11, 0.65]: -11158.852689423227 [9.75, -1.51, 1083.61, 0.1, -3.11, 0.27, 0.25]: -8993.678119508962 [6.77, -0.99, 2213.42, 0.03, -3.21, 0.06, 0.31]: -13885.689933092046 [9.0, -2.73, 2310.76, 0.06, -1.12, 0.76, 0.69]: -12218.209174197544
rand_pars = [9.5, 0.0, 500, 0.1, -2, 0.5, 0.4]
model.lnpost(rand_pars)
-4812.512702328114
import pickle
pickle.dump(model, open('testmodel.pkl', 'wb'))
model = pickle.load(open('testmodel.pkl', 'rb'))
model.lnpost(pars)
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) <ipython-input-11-376fe96b2274> in <module>() ----> 1 model.lnpost(pars) /Users/tdm/repositories/isochrones/isochrones/starmodel.py in lnpost(self, p, **kwargs) 487 if not np.isfinite(lnpr): 488 return lnpr --> 489 return lnpr + self.lnlike(p, **kwargs) 490 491 def lnlike(self, p, **kwargs): /Users/tdm/repositories/isochrones/isochrones/cluster.py in lnlike(self, p) 147 148 # Compute model mags at each eep --> 149 model_mags = {b : self.ic.mag[b](eeps, age, feh, distance, AV) for b in self.bands} 150 151 # Compute the log-likelihood of the (non-mag) props /Users/tdm/repositories/isochrones/isochrones/cluster.py in <dictcomp>(.0) 147 148 # Compute model mags at each eep --> 149 model_mags = {b : self.ic.mag[b](eeps, age, feh, distance, AV) for b in self.bands} 150 151 # Compute the log-likelihood of the (non-mag) props /Users/tdm/repositories/isochrones/isochrones/isochrone.py in fn(mass, age, feh, distance, AV, x_ext, ext_table) 245 A = AV*ext(LAMBDA_EFF[band]) 246 dm = 5*np.log10(distance) - 5 --> 247 return self._mag[band](mass, age, feh) + dm + A 248 return fn 249 AttributeError: 'MIST_Isochrone' object has no attribute '_mag'
model.fit(overwrite=True)
/Users/tdm/repositories/isochrones/isochrones/cluster.py:119: RuntimeWarning: divide by zero encountered in log val = np.log(self._priors[prop](eval(prop))) ERROR: Interrupt received: Terminating
! say "sampling complete!"
%matplotlib inline
import corner
names = ['age', 'feh', 'distance', 'AV', 'gamma']
fig = corner.corner(model.samples[names], names=names,
truths=[age, feh, distance, AV, gamma])
WARNING:root:Too few points to create valid contours WARNING:root:Too few points to create valid contours
bands = 'rJK'
def band_pairs(bands):
return [(bands[i], bands[-1]) for i in range(len(bands)-1)]
band_pairs(bands)
[('r', 'K'), ('J', 'K')]
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)
%%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)
DynamicMap cannot be displayed without explicit indexing as 'AV1', 'AV2' dimension(s) are unbounded. Set dimensions bounds with the DynamicMap redim.range or redim.values methods.
:DynamicMap [age1,feh1,dist1,AV1,age2,feh2,dist2,AV2]
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)
age2, feh2, dist2, AV2, gamma2 = other_pars
iso_compare(mist, bands='rJK', age1=age, feh1=feh, dist1=distance, AV1=AV,
age2=age2, feh2=feh2, dist2=dist2, AV2=AV2)
means = model.samples.mean()
means
age 9.493336 feh -0.226858 distance 49877.166757 AV 0.002651 gamma -1.358348 lnprob -54.972238 dtype: float64
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)
(-11.165645257906878, -52.74261548012587, -63.90826073803275)
model.lnprior(pars), model.lnlike(pars), model.lnpost(pars)
(-16.641241209353648, -186.79573788050948, -203.43697908986312)
%debug model.lnlike(pars)
NOTE: Enter 'c' at the ipdb> prompt to continue execution. > <string>(1)<module>() ipdb> s --Call-- > /Users/tdm/repositories/isochrones/isochrones/cluster.py(76)lnlike() 74 return lnp 75 ---> 76 def lnlike(self, p): 77 age = p[0] 78 feh = p[1] ipdb> n > /Users/tdm/repositories/isochrones/isochrones/cluster.py(77)lnlike() 75 76 def lnlike(self, p): ---> 77 age = p[0] 78 feh = p[1] 79 distance = p[2] ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(78)lnlike() 76 def lnlike(self, p): 77 age = p[0] ---> 78 feh = p[1] 79 distance = p[2] 80 AV = p[3] ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(79)lnlike() 77 age = p[0] 78 feh = p[1] ---> 79 distance = p[2] 80 AV = p[3] 81 gamma = p[4] ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(80)lnlike() 78 feh = p[1] 79 distance = p[2] ---> 80 AV = p[3] 81 gamma = p[4] 82 ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(81)lnlike() 79 distance = p[2] 80 AV = p[3] ---> 81 gamma = p[4] 82 83 lnlike_tot = 0 ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(83)lnlike() 81 gamma = p[4] 82 ---> 83 lnlike_tot = 0 84 85 mineep, maxeep = self.bounds('eep') ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(85)lnlike() 83 lnlike_tot = 0 84 ---> 85 mineep, maxeep = self.bounds('eep') 86 eeps = np.arange(mineep, maxeep + 1) 87 ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(86)lnlike() 84 85 mineep, maxeep = self.bounds('eep') ---> 86 eeps = np.arange(mineep, maxeep + 1) 87 88 # Compute log-likelihood of each mass under power-law distribution ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(90)lnlike() 88 # Compute log-likelihood of each mass under power-law distribution 89 # Also use this opportunity to find the valid range of EEP ---> 90 mass_fn = PowerLawPrior(gamma, bounds=(self.ic.minmass, self.ic.maxmass)) 91 model_masses = self.ic.initial_mass(eeps, age, feh) 92 ok = np.isfinite(model_masses) ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(91)lnlike() 89 # Also use this opportunity to find the valid range of EEP 90 mass_fn = PowerLawPrior(gamma, bounds=(self.ic.minmass, self.ic.maxmass)) ---> 91 model_masses = self.ic.initial_mass(eeps, age, feh) 92 ok = np.isfinite(model_masses) 93 ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(92)lnlike() 90 mass_fn = PowerLawPrior(gamma, bounds=(self.ic.minmass, self.ic.maxmass)) 91 model_masses = self.ic.initial_mass(eeps, age, feh) ---> 92 ok = np.isfinite(model_masses) 93 94 model_masses = model_masses[ok] ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(94)lnlike() 92 ok = np.isfinite(model_masses) 93 ---> 94 model_masses = model_masses[ok] 95 eeps = eeps[ok] 96 ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(95)lnlike() 93 94 model_masses = model_masses[ok] ---> 95 eeps = eeps[ok] 96 97 lnlike_mass = np.log(mass_fn.pdf(model_masses)) ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(97)lnlike() 95 eeps = eeps[ok] 96 ---> 97 lnlike_mass = np.log(mass_fn.pdf(model_masses)) 98 99 # Compute log-likelihood of observed photometry ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(100)lnlike() 98 99 # Compute log-likelihood of observed photometry --> 100 model_mags = {b : self.ic.mag[b](eeps, age, feh, distance, AV) for b in self.bands} 101 102 lnlike_phot = 0 ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(102)lnlike() 100 model_mags = {b : self.ic.mag[b](eeps, age, feh, distance, AV) for b in self.bands} 101 --> 102 lnlike_phot = 0 103 for b in self.bands: 104 vals = self.stars[b].values ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(103)lnlike() 101 102 lnlike_phot = 0 --> 103 for b in self.bands: 104 vals = self.stars[b].values 105 uncs = self.stars[b + '_unc'].values ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(104)lnlike() 102 lnlike_phot = 0 103 for b in self.bands: --> 104 vals = self.stars[b].values 105 uncs = self.stars[b + '_unc'].values 106 ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(105)lnlike() 103 for b in self.bands: 104 vals = self.stars[b].values --> 105 uncs = self.stars[b + '_unc'].values 106 107 lnlike_phot += -0.5 * (vals - model_mags[b][:, None])**2 / uncs**2 ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(107)lnlike() 105 uncs = self.stars[b + '_unc'].values 106 --> 107 lnlike_phot += -0.5 * (vals - model_mags[b][:, None])**2 / uncs**2 108 109 integrand = np.exp(lnlike_mass[:, None] + lnlike_phot) ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(103)lnlike() 101 102 lnlike_phot = 0 --> 103 for b in self.bands: 104 vals = self.stars[b].values 105 uncs = self.stars[b + '_unc'].values ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(104)lnlike() 102 lnlike_phot = 0 103 for b in self.bands: --> 104 vals = self.stars[b].values 105 uncs = self.stars[b + '_unc'].values 106 ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(105)lnlike() 103 for b in self.bands: 104 vals = self.stars[b].values --> 105 uncs = self.stars[b + '_unc'].values 106 107 lnlike_phot += -0.5 * (vals - model_mags[b][:, None])**2 / uncs**2 ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(107)lnlike() 105 uncs = self.stars[b + '_unc'].values 106 --> 107 lnlike_phot += -0.5 * (vals - model_mags[b][:, None])**2 / uncs**2 108 109 integrand = np.exp(lnlike_mass[:, None] + lnlike_phot) ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(103)lnlike() 101 102 lnlike_phot = 0 --> 103 for b in self.bands: 104 vals = self.stars[b].values 105 uncs = self.stars[b + '_unc'].values ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(104)lnlike() 102 lnlike_phot = 0 103 for b in self.bands: --> 104 vals = self.stars[b].values 105 uncs = self.stars[b + '_unc'].values 106 ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(105)lnlike() 103 for b in self.bands: 104 vals = self.stars[b].values --> 105 uncs = self.stars[b + '_unc'].values 106 107 lnlike_phot += -0.5 * (vals - model_mags[b][:, None])**2 / uncs**2 ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(107)lnlike() 105 uncs = self.stars[b + '_unc'].values 106 --> 107 lnlike_phot += -0.5 * (vals - model_mags[b][:, None])**2 / uncs**2 108 109 integrand = np.exp(lnlike_mass[:, None] + lnlike_phot) ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(103)lnlike() 101 102 lnlike_phot = 0 --> 103 for b in self.bands: 104 vals = self.stars[b].values 105 uncs = self.stars[b + '_unc'].values ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(109)lnlike() 107 lnlike_phot += -0.5 * (vals - model_mags[b][:, None])**2 / uncs**2 108 --> 109 integrand = np.exp(lnlike_mass[:, None] + lnlike_phot) 110 111 like_tot = np.trapz(integrand, axis=0) ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(111)lnlike() 109 integrand = np.exp(lnlike_mass[:, None] + lnlike_phot) 110 --> 111 like_tot = np.trapz(integrand, axis=0) 112 113 # Don't log a zero. ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(114)lnlike() 112 113 # Don't log a zero. --> 114 ok = (like_tot != 0) 115 if not ok.any(): 116 return -np.inf ipdb> > /Users/tdm/repositories/isochrones/isochrones/cluster.py(115)lnlike() 113 # Don't log a zero. 114 ok = (like_tot != 0) --> 115 if not ok.any(): 116 return -np.inf 117 lnlike = np.log(like_tot[ok]).sum() ipdb> like_tot array([3.21015311e-03, 2.29907191e-03, 1.48631766e-02, 3.31035104e-04, 2.93435232e-03, 1.34640436e-03, 4.68624678e-03, 8.97556563e-04, 8.93346158e-04, 7.30555379e-03, 1.06285935e-03, 8.62753872e-03, 1.92270572e-03, 3.11094086e-03, 1.66123822e-03, 8.67627177e-03, 8.62511112e-04, 7.67518574e-04, 5.79266336e-03, 1.48829421e-03, 6.67204447e-04, 3.68704256e-03, 1.63487285e-03, 1.86675657e-03, 7.03456281e-05, 1.66416498e-03, 7.50282216e-04, 2.06727392e-03, 9.40632898e-03, 3.29603662e-03]) ipdb> q