In [1]:
%load_ext autoreload
%autoreload 2
In [1]:
%matplotlib
import matplotlib.pyplot as plt
import numpy as np
from IPython.html import widgets # Widget definitions
from IPython.display import display, clear_output, HTML # Used to display widgets in the notebook
from IPython.html.widgets import interact, interactive, fixed
In [2]:
import Starfish
import Starfish.grid_tools
from Starfish.samplers import StateSampler
from Starfish.spectrum import DataSpectrum, Mask, ChebyshevSpectrum
from Starfish.emulator import Emulator
import Starfish.constants as C
from Starfish.covariance import get_dense_C, make_k_func, make_k_func_region
from Starfish.model import ThetaParam, PhiParam
from scipy.special import j1
from scipy.interpolate import InterpolatedUnivariateSpline
import numpy as np
#myDataSpectrum = DataSpectrum.open("../data/WASP14/WASP14-2010-03-29.hdf5", orders=np.array([22]))
myDataSpectrum = DataSpectrum.open("../data/WASP14/WASP14-2009-06-14.hdf5", orders=np.array([22]))
#mask = np.load("WASP14_24.mask.npy")
#myDataSpectrum.add_mask(np.atleast_2d(mask))
#myDataSpectrum = DataSpectrum.open("../data/LkCa15/LkCa15-2013-10-12.hdf5", orders=np.array([32]))
# myDataSpectrum = DataSpectrum.open("../data/Gl51/Gl51.hdf5")
myInstrument = TRES()
# myInstrument = SPEX()
myHDF5Interface = HDF5Interface("../libraries/PHOENIX_TRES_F.hdf5")
#myHDF5Interface = HDF5Interface("../libraries/Kurucz_TRES.hdf5")
#myHDF5Interface = HDF5Interface("../libraries/Kurucz_master.hdf5")
In [13]:
#Load a model using the JSON file
myModel = Model.from_json("WASP2_21_22_23_model_final.json", myDataSpectrum, myInstrument, myHDF5Interface)
#myModel.evaluate()
In [3]:
#Or, instantiate and tweak parameters
myModel = Model(myDataSpectrum, myInstrument, myHDF5Interface, stellar_tuple=("temp", "logg", "Z", "vsini", "vz", "logOmega"),
cheb_tuple=("logc0","c1", "c2", "c3"), cov_tuple=("sigAmp", "logAmp", "l"), region_tuple=("loga", "mu", "sigma"),
max_v=20, debug=True)
In [4]:
#myOrderModel.update_Cov({"sigAmp": 1, "logAmp":-14 , "l":10.})
#First pass at actually plotting the model
#params = {"temp":6000, "logg":4.0, "Z":-0.3, "vsini":4, "vz":-5, "logOmega":-19.6}
#params = {"temp":6315, "logg":3.92, "Z":-0.3, "vsini":4, "vz":-5, "logOmega":-12.72}
#params = {'temp': 6313.276253398848, 'vz': -4.5886794506100905, 'logOmega': -19.659,
# 'Z': -0.45755240763502216, 'vsini': 5.1871494360273905, 'logg': 3.9089865093046496}
params = {'vz': -4.8456827988388413, 'logOmega': -19.658942686778765, 'vsini': 5.0328046854058668,
'temp': 6220.0106396310739, 'logg': 4.29, 'Z': -0.35001081006631218}
#params = {"temp":6315, "logg":3.92, "Z":-0.3, "vsini":4, "vz":-4.4, "logOmega":-12.96}
#params = {"temp":4500, "logg":4.0, "Z":-0.1, "vsini":10, "vz":15, "logOmega":-19.}
# params = {"temp":3000, "logg":4.0, "Z":-0.1, "vsini":10, "vz":0, "logOmega":-19.}
myModel.update_Model(params)
print(myModel.evaluate())
#myOrderModel.update_Cheb({"c1":0.0, "c2":0.0, "c3":0.0})
#cov = myModel.OrderModels[0].get_Cov()
In [5]:
#params = {'logOmega': -19.69045147013749, 'temp': 6219.4026633831281, 'vz': -4.864716315534733, 'logg': 4.29,
# 'Z': -0.25012131204369881, 'vsini': 4.994423918348339}
#params = {'logOmega': -19.659045147013749, 'temp': 6219.4026633831281, 'vz': -4.864716315534733, 'logg': 4.29,
# 'Z': -0.3509729219906163, 'vsini': 4.994423918348339}
params = {'vz': -4.8456827988388413, 'logOmega': -19.658942686778765, 'vsini': 5.0328046854058668, 'temp': 6220.0106396310739,
'logg': 4.29, 'Z': -0.35001081006631218}
myModel.update_Model(params)
print(myModel.evaluate())
#Strangely, Z is totally the value that changes this. Why? This is super-precise inference on Z, to the 4th digit. WTF?
In [63]:
myOrderModel = myModel.OrderModels[0]
#myOrderModel.update_Cheb({'c3': 0.0022745611258502803, 'c2': -0.024660521505542141,
# 'logc0': -0.020399928159079816, 'c1': -0.05817214397082078})
myOrderModel.update_Cheb({'c3': 0.00, 'c2': 0.0,
'logc0': 0.0, 'c1': -0.0})
print(myOrderModel.evaluate())
#myOrderModel = myModel.OrderModels[1]
model_flux = myOrderModel.get_spectrum()
In [10]:
spec = myModel.get_data()
wl = spec.wls[0]
fl = spec.fls[0]
mask = spec.masks[0]
#wl = spec.wls[1]
#fl = spec.fls[1]
#mask = spec.masks[1]
fig, ax = plt.subplots(nrows=2, figsize=(11,8), sharex=True)
ax[0].plot(wl, fl, "b")
l_model, = ax[0].plot(wl, model_flux, "r")
ax[0].set_ylabel("Data and Model")
l_resid, = ax[1].plot(wl, fl - model_flux, "b")
ax[1].set_xlabel(r"$\lambda$\AA")
ax[1].set_ylabel("Residuals")
#good = [5142.8, 5150.7, 5151.8, 5169.1, 5170.7, 5172.3, 5183.2, 5183.8, 5184.2, 5187.8, 5188.8, 5197.4, 5198.6, 5217.7, 5235.3]
#bad = [5145., 5167.0, 5171.5, 5178.7, 5194.9, 5202.3, 5226.5, 5232.5, 5234.6]
#Plot the axvline on the lines
#for line in good:
# ax[1].axvline(line, color="0.5")
#for line in bad:
# ax[1].axvline(line, color="r")
#cov = myModel.OrderModels[0].get_Cov().todense()
#fig2 = plt.figure()
#ax2 = fig2.add_subplot(111)
#im = ax2.imshow(cov, origin='upper', interpolation='none')
plt.show()
In [9]:
def update_model_plot(**kwargs):
'''Take the kwargs, update the model and residuals'''
#Update the model spectrum
myModel.update_Model(kwargs)
model_flux = myOrderModel.get_spectrum()
l_model.set_ydata(model_flux)
#Update the residuals
residuals = fl - model_flux
l_resid.set_ydata(residuals)
#Find ymax and ymin and rescale
ax[0].set_ylim(np.min([fl, model_flux]), np.max([fl, model_flux]))
ax[1].set_ylim(np.min(fl - model_flux), np.max(fl - model_flux))
#Redraw the plot
fig.canvas.draw_idle()
#Calculate and print the lnprob
print(myOrderModel.evaluate())
In [11]:
np.save("WASPfl.npy", myOrderModel.get_spectrum())
np.save("WASP_resid.npy", fl - model_flux)
In [8]:
def update_Cheb_plot(**kw):
'''Take the kwargs, update the model and residuals'''
#Update the Chebyshev polynomial
myOrderModel.update_Cheb(kw)
model_flux = myOrderModel.get_spectrum()
l_model.set_ydata(model_flux)
#Update the residuals
l_resid.set_ydata(fl - model_flux)
#Find ymax and ymin and rescale
ax[0].set_ylim(np.min([fl, model_flux]), np.max([fl, model_flux]))
ax[1].set_ylim(np.min(fl - model_flux), np.max(fl - model_flux))
#Redraw the plot
fig.canvas.draw_idle()
#Calculate and print the lnprob
#print(myModel.evaluate())
In [8]:
def update_Cov_plot(**kwargs):
'''Take the kwargs, update the model and residuals'''
#Update the covariance matrix
myModel.OrderModels[0].update_Cov(kwargs)
cov = myModel.OrderModels[0].get_Cov().todense()
#Replot the covariance matrix
im.set_array(cov)
#Redraw the plot
fig2.canvas.draw_idle()
#Calculate and print the lnprob
print(myModel.evaluate())
In [11]:
i = interact(update_model_plot,
temp=(5500,6500, 10),
logg=(2.5,5.0, 0.1),
Z=(-1, 0.5, 0.1),
#alpha=(0.0, 0.4, 0.05),
vsini=(3, 8., 0.5),
vz=(-10,0, 0.1),
#Av=(0.0,1.0, 0.05),
logOmega=(-13.,-12., 0.01),
#logOmega=(-20.4,-19.0, 0.01),
)
In [10]:
i = interact(update_Cheb_plot,
c1=(-0.2, 0.2, 0.01),
c2=(-0.2, 0.2, 0.01),
c3=(-0.2, 0.2, 0.01),
)
In [14]:
i = interact(update_Cov_plot,
sigAmp=(.5,1.5, 0.1),
logAmp=(-15,-13, 0.2),
l=(1, 50, 1),
)
Good starting guesses for WASP-14: temp: 6100 logg: 4.0 Z: -0.5 vsini: 6 vz: 13.7 log_Omega: -19.7 alpha: 0.2