In [1]:
%load_ext autoreload
%autoreload 2
In [2]:
#Import all of the relevant python modules
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter as FSF
from matplotlib.ticker import MaxNLocator
from matplotlib.ticker import MultipleLocator
In [3]:
from Starfish.model import Model
from Starfish.spectrum import DataSpectrum
from Starfish.grid_tools import TRES, HDF5Interface
from Starfish import utils
from Starfish.utils import saveall
import scipy.sparse as sp
import numpy as np
myDataSpectrum = DataSpectrum.open("../../data/WASP14/WASP14-2009-06-14.hdf5", orders=np.array([22]))
myInstrument = TRES()
myHDF5Interface = HDF5Interface("../../libraries/PHOENIX_TRES_F.hdf5")
In [3]:
#Load a model using the JSON file
#myModel = Model.from_json("../WASP14/model_final.json", myDataSpectrum, myInstrument, myHDF5Interface)
myModel = Model.from_json("../WASP14_22_model_final_region.json", myDataSpectrum, myInstrument, myHDF5Interface)
myOrderModel = myModel.OrderModels[0]
spec = myModel.get_data()
wl = spec.wls[0]
fl = spec.fls[0]
sigma = spec.sigmas[0]
model_fl = myOrderModel.get_spectrum()
residuals = fl - model_fl
cheb = myOrderModel.get_Cheb()
In [6]:
sh0 = 5201.4
sh1 = 5205.1
In [1]:
fig, ax = plt.subplots(nrows=2, figsize=(7,2.05), sharex=True)
ax[0].plot(wl, fl*1e13, "w", lw=1.8)
ax[0].plot(wl, model_fl*1e13, "w", lw=1.8)
l0, = ax[0].plot(wl, fl*1e13, "b")
l1, = ax[0].plot(wl, model_fl*1e13, "r")
#ax[0].legend([l0, l1], ["data", "model"], loc="lower right", ncol=2, prop={'size':10})
ax[0].yaxis.set_major_locator(MaxNLocator(nbins=7, prune='lower'))
ax[0].fill_betweenx(np.array([0.1, 2.4]), sh0, sh1, color="0.5", alpha=0.5)
ax[0].set_ylim(0.1, 2.3)
ax[0].annotate("data", (0.86, 0.1), xycoords="axes fraction", ha="right", color="b", size=9)
ax[0].annotate("model", (0.97, 0.1), xycoords="axes fraction", ha="right", color="r", size=9)
ax[1].fill_betweenx(np.array([-0.9, 0.9]), sh0, sh1, color="0.5", alpha=0.5)
ax[1].plot(wl, residuals*1e13, "w", lw=1.8)
l2, = ax[1].plot(wl, residuals*1e13, "k")
ax[1].xaxis.set_major_formatter(FSF("%.0f"))
ax[1].set_xlabel(r"$\lambda$ [\AA]")
for a in ax:
a.yaxis.set_major_locator(MaxNLocator(nbins=5, prune='both'))
ax[1].set_xlim(5190, 5236)
ax[1].set_ylim(-0.45, 0.45)
#ax[1].legend([l2], ["residuals"], loc="upper right", prop={'size':10})
ax[1].annotate("residuals", (0.97, 0.8), xycoords="axes fraction", ha="right", color="k", size=9)
fig.subplots_adjust(hspace=0, left=0.12, right=0.88, bottom=0.19, top=0.94)
fig.text(0.035, 0.7, r"$f_\lambda\, \times 10^{-13}$" + "\n"
+ r"$[\textrm{erg}\;\textrm{cm}^{-2}\;\textrm{s}^{-1}\;\textrm{\AA}^{-1}]$", rotation="vertical", ha="center")
saveall(fig, "../../plots/matrix_residual")
plt.show()
In [8]:
wl = myDataSpectrum.wls[0]
sigma = myDataSpectrum.sigmas[0]
In [9]:
ind = (wl > sh0) & (wl < sh1)
w = wl[ind]
s = sigma[ind]
#r = residuals[ind]
In [11]:
#These are using some of the best-fit parameters from
#/n/home07/iczekala/StellarSpectra/output/WASP14/PHOENIX/22/region/4sig/2014-08-13/run60/22
mat_poisson = utils.Poisson_matrix(w, sigma=s)
mat_global = utils.k_global_matrix(w, a=10**(-14.34), l=10.26)
# -13.77186598, 5202.25227836, 7.45325824
mat_local = utils.k_local_matrix(w, a=10**(-13.77), mu=5202.25, sigma=7.4532)
mat_combined = mat_poisson + mat_global + mat_local
In [28]:
more_draws = utils.random_draws(mat_combined, num=200)
In [15]:
draws = utils.random_draws(mat_combined, num=100)
np.save("combined_draws.npy", draws)
In [29]:
draws = np.load("combined_draws.npy")
In [16]:
utils.visualize_draws(draws)
In [30]:
#Colorbrewer bands
s3 = '#fee6ce'
s2 = '#fdae6b'
s1 = '#e6550d'
In [113]:
def miniplot(ax, cov, spaces):
'''
Given an axes, plot the scenario of a bunch of random draws on top of it
'''
draws = utils.random_draws(cov, num=200)
#plot the first three staggered
#spaces = sep = 6*np.std(draws[0])
for draw, space in zip(draws[:3], spaces):
ax.axhline(space*1e13, ls=":", color="0.3")
ax.plot(w, (draw + space)*1e13, "DarkGreen", ls="steps-mid")
#Jump again and now plot the envelope of the draws
min_spec, max_spec = utils.std_envelope(draws)
ax.fill_between(w, 3*min_spec*1e13, 3*max_spec*1e13, zorder=0, color=s3)
ax.fill_between(w, 2*min_spec*1e13, 2*max_spec*1e13, zorder=0, color=s2)
ax.fill_between(w, min_spec*1e13, max_spec*1e13, zorder=0, color=s1)
ax.plot(w, r*1e13, "w", lw=1.8, zorder=1, ls="steps-mid")
ax.plot(w, r*1e13, "k", zorder=1, ls="steps-mid")
In [114]:
def complot(ax, draws, selected, spaces):
'''
Given selected draws, copy the plot
'''
for draw, space in zip(selected, spaces):
ax.axhline(space*1e13, ls=":", color="0.3")
ax.plot(w, (draw + space)*1e13, "DarkGreen", ls="steps-mid")
#Jump again and now plot the envelope of the draws
min_spec, max_spec = utils.std_envelope(draws)
ax.fill_between(w, 3*min_spec*1e13, 3*max_spec*1e13, zorder=0, color=s3)
ax.fill_between(w, 2*min_spec*1e13, 2*max_spec*1e13, zorder=0, color=s2)
ax.fill_between(w, min_spec*1e13, max_spec*1e13, zorder=0, color=s1)
ax.plot(w, r*1e13, "w", lw=1.8, zorder=1, ls="steps-mid")
ax.plot(w, r*1e13, "k", zorder=1, ls="steps-mid")
In [115]:
fig, ax = plt.subplots(nrows=3, figsize=(3.2, 6.), sharex=True)
spaces = 1.1 * np.array([3.97828569937e-14, 7.95657139874e-14, 1.19348570981e-13]) #determined from
#spaces = sep = 6*np.std(draws[0]) on a draw from the combined matrix
miniplot(ax[0], mat_poisson, spaces)
miniplot(ax[1], mat_global + mat_poisson, spaces)
complot(ax[2], draws=more_draws, selected=draws[np.array([17, 6, 16])], spaces=spaces)
l0, = ax[0].plot(np.ones(10), "0.4")
l1, = ax[0].plot(np.ones(10), "k")
#legend1 = ax[0].legend([l0], ["random draw"], loc="upper right", prop={'size':10})
#ax[0].legend([l1], ["residuals"], loc="lower right", prop={'size':10})
#ax[0].add_artist(legend1)
ax[0].annotate("random draw", (0.97, 0.9), xycoords="axes fraction", ha="right", color="DarkGreen", size=9)
ax[0].annotate("residuals", (0.97, 0.09), xycoords="axes fraction", ha="right", color="k", size=9)
fig.text(0.03, 0.9, r"$f_\lambda$ + offset $\times 10^{-13}\,[\textrm{erg}\;\textrm{cm}^{-2}\;\textrm{s}^{-1}\;\textrm{\AA}^{-1}]$",
rotation="vertical")
ax[-1].set_xlim(sh0 + 0.1, sh1 - 0.1)
ax[-1].xaxis.set_major_formatter(FSF("%.0f"))
ax[-1].xaxis.set_major_locator(MultipleLocator(1.))
ax[-1].set_xlabel(r"$\lambda$ [\AA]")
for a in ax:
a.yaxis.set_major_locator(MaxNLocator(nbins=7, prune='upper'))
a.set_ylim(-.55, 1.75)
for a in ax[:2]:
a.yaxis.set_ticklabels([])
fig.subplots_adjust(hspace=0.05, top=0.98, right=0.94, left=0.13, bottom=0.08)
saveall(fig, "../../plots/white_to_local")
plt.show()
In [43]:
np.log10(np.max(mat_poisson + mat_global + mat_local))
Out[43]:
In [41]:
plt.hist(mat_poisson.flatten(), bins=100)
plt.show()
In [ ]:
np.log10()
In [13]:
norm = matplotlib.colors.Normalize(vmin=-31, vmax=-27.5, clip=True)
In [14]:
def plot_matrix(wl, cov, title=False, axl=False):
'''
Construct the plot!
'''
fig = plt.figure(figsize=(2.95,2.35))
ax = fig.add_subplot(111)
left, right = wl[0], wl[-1]
cov += 1e-99
img = ax.imshow(np.log10(cov), interpolation="none", origin="upper", extent=[left, right, right, left], cmap=plt.cm.Blues, norm=norm)
ax.xaxis.set_major_formatter(FSF("%.0f"))
ax.xaxis.set_major_locator(MultipleLocator(1.))
#labels = ax.get_xticklabels()
#for label in labels:
# label.set_rotation(40)
ax.yaxis.set_major_formatter(FSF("%.0f"))
ax.yaxis.set_major_locator(MultipleLocator(1.))
if axl:
ax.set_xlabel(r"$\lambda$ [\AA]")
ax.set_ylabel(r"$\lambda$ [\AA]")
else:
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
if title:
cax = fig.add_axes([0.82, 0.22, 0.03, 0.65])
cb = fig.colorbar(img, cax=cax)
#ticks = np.linspace(0, np.max(cov), num=6)
#cb.set_ticks(ticks)
cb.set_ticks(MaxNLocator(nbins=5))
fig.text(0.25, 0.9, r"$\mathsf{C}_{ij}\,\log_{10}\,[\textrm{erg}^2\;\textrm{cm}^{-4}\;\textrm{s}^{-2}\;\textrm{\AA}^{-2}]$")
fig.subplots_adjust(left=0.2, bottom=0.1, top=0.98, right=0.79)
return fig,ax
#saveall(fig, fname)
#plt.show()
In [15]:
fig, ax = plot_matrix(w, mat_poisson, title=True, axl=False)
ax.annotate(r"$\delta_{ij} \sigma_i$", (5204.95, 5201.95), ha="right", size=9)
saveall(fig, "../../plots/matrix_poisson")
plt.show()
In [16]:
fig, ax = plot_matrix(w, mat_poisson + mat_global, axl=False)
ax.annotate(r"$\mathcal{K}^{\rm G} + \delta_{ij} \sigma_i$", (5204.95, 5201.95), ha="right", size=9)
saveall(fig, "../../plots/matrix_global")
plt.show()
In [17]:
fig, ax = plot_matrix(w, mat_poisson + mat_global + mat_local, axl=True)
ax.annotate(r"$\mathcal{K}^{\rm L} + \mathcal{K}^{\rm G} + \delta_{ij} \sigma_i$", (5204.95, 5201.95), ha="right", size=9)
saveall(fig, "../../plots/matrix_local")
plt.show()
In [ ]: