In [99]:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.ensemble import ExtraTreesRegressor
from time import time
from scipy import optimize
import os
import ipywidgets
from ipywidgets import interact, interactive, fixed
In [2]:
simulations_filename = os.path.join('..', 'forward', 'simulations.dat')
data = pd.read_csv(simulations_filename, sep='\t')
In [3]:
data.head()
Out[3]:
In [4]:
indices = np.ndarray.flatten(np.array([range(64) for i in range(int(len(data)/64))]))
data['index'] = indices
data.head()
Out[4]:
In [5]:
X_names = ['M', 'Y', 'Z', 'alpha', 'overshoot', 'diffusion', 'index']
X = data.loc[:,X_names]
X.head()
Out[5]:
In [6]:
y_names = ['age', 'X_c', 'radius', 'L', 'log_g', 'Teff', 'Fe/H', 'Dnu0', 'dnu02']
y = data.loc[:,y_names]
y.head()
Out[6]:
In [142]:
forest = ExtraTreesRegressor(n_estimators=64, n_jobs=1, oob_score=True, bootstrap=True)
start = time()
forest.fit(X, y)
end = time()
print(forest.oob_score_, end-start)
In [24]:
np.vstack((X.min(), X.max())).T
Out[24]:
In [146]:
inputs = [1, 0.2745, 0.02, 1.9139, 0.0518, 1, 29.3161]
res = forest.predict(np.array(inputs).reshape(1, -1))[0]
print(res[y_names.index('Teff')])
print(res[y_names.index('L')])
print(res[y_names.index('age')])
print(res[y_names.index('radius')])
In [141]:
def Teff_L_age(x):
Y, Z, alpha_MLT, alpha_ov, index = x
#Y, alpha_MLT = x
inputs = [1, Y, Z, alpha_MLT, alpha_ov, 1, index]
#inputs = [1, Y, 0.02, alpha_MLT, 0.05, 1, 29]
res = forest.predict(np.array(inputs).reshape(1, -1))[0]
return (res[y_names.index('Teff')] - 5777)**2 / 100 + \
(res[y_names.index('L')] - 1)**2 + \
(res[y_names.index('age')] - 4.57)**2/0.3 + \
(res[y_names.index('radius')] - 1)**2
result = optimize.minimize(Teff_L_age, x0=[0.275, 0.02, 1.9, 0.05, 29], method='Nelder-Mead')
In [143]:
print(list(map(lambda x: "{0:.4f}".format(x), result.x)))
In [144]:
def plot_HR(M, Y, Z, alpha_MLT, alpha_ov, diffusion): #, index):
plt.ion()
#x=('apples','oranges')
Teff = []
L = []
for ii in range(63):
res = forest.predict(np.array([M, Y, Z, alpha_MLT, alpha_ov, diffusion, ii]).reshape(1, -1))[0]
Teff += [res[y_names.index('Teff')]]
L += [res[y_names.index('L')]]
plt.plot(Teff, L, 'k-')
plt.gca().invert_xaxis()
plt.text(5777, 1, r'$\odot$')
plt.xlabel(r'Effective Temperature $T_{\mathrm{eff}}/K$')
plt.ylabel(r'Luminosity $L/L_\odot$')
return
def interactive():
M_slider = ipywidgets.FloatSlider(min=0.7, max=1.6, step=0.01, value=1)
Y_slider = ipywidgets.FloatSlider(min=0.22, max=0.34, step=0.001, value=0.2745)
Z_slider = ipywidgets.FloatSlider(min=10e-5, max=0.1, step=0.0001, value=0.02)
a_slider = ipywidgets.FloatSlider(min=1.5, max=2.5, step=0.01, value=1.9139)
o_slider = ipywidgets.FloatSlider(min=0, max=0.1, step=0.01, value=0.0518)
D_slider = ipywidgets.FloatSlider(min=0, max=100, step=0.1, value=1)
interact(plot_HR, M=M_slider, Y=Y_slider, Z=Z_slider, alpha_MLT=a_slider, alpha_ov=o_slider,
diffusion=D_slider)#, index=i_slider)
In [145]:
%matplotlib inline
interactive()