In [1]:

    

from Starfish.emulator import PCAGrid
from Starfish.grid_tools import HDF5Interface
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.ticker import FormatStrFormatter as FSF
from matplotlib.ticker import MaxNLocator
from matplotlib.ticker import MultipleLocator
from Starfish.utils import saveall


    

In [2]:

    

pca = PCAGrid.open("../../libraries/PHOENIX_SPEX_M_PCA.hdf5")
interface = HDF5Interface("../../libraries/PHOENIX_SPEX_M.hdf5")
wl = pca.wl


    

In [3]:

    

ind = (wl[0] <= interface.wl) & (interface.wl <= wl[-1])


    

In [4]:

    

rfluxes = np.array([flux[ind] for flux in interface.fluxes])
fluxes = np.array([flux[ind] for flux in interface.fluxes])

#Normalize all of the fluxes to an average value of 1
#In order to remove interesting correlations
flux_avg = np.average(fluxes, axis=1)[np.newaxis].T

fluxes = fluxes/flux_avg
rfluxes = rfluxes/flux_avg

#Subtract the mean from all of the fluxes.
flux_mean = np.average(fluxes, axis=0)
fluxes -= flux_mean

#"Whiten" each spectrum such that the variance for each wavelength is 1
flux_std = np.std(fluxes, axis=0)
fluxes /= flux_std


    

In [5]:

    

interface.list_grid_points[187]


    

    
Out[5]:





{'logg': 5.0, 'Z': 0.0, 'alpha': 0.0, 'temp': 3000}

In [6]:

    

pca.gparams[37]


    

    
Out[6]:





array([ 3000.,     5.,     0.])

In [70]:

    

#h_spec = interface.load_flux({"temp":3000, "logg":5.0, "Z":0.0, "alpha":0.0})[ind]


    

In [7]:

    

h_spec = rfluxes[187]


    

In [8]:

    

std_spec = fluxes[187]


    

In [9]:

    

raw_spec = fluxes[187] * flux_std + flux_mean


    

In [10]:

    

ww = pca.pcomps.dot(std_spec)


    

In [11]:

    

ww


    

    
Out[11]:





array([-27.89909803, -10.62579203,  -0.66197484,  -3.63912344,  -2.35353456])

In [13]:

    

recon_spec = pca.reconstruct(ww)


    

In [35]:

    

fig, ax = plt.subplots(nrows=2, figsize=(7,4.3), sharex=True)

wlm = wl * 1e-4 #convert to microns

w_list = [r"w_{}".format(j) for j in range(1,6)]
xi_list = [r"$\xi_{}$".format(j) for j in range(1,6)]

ax[0].plot(wlm, pca.flux_mean + 0.9, "0.5")
ax[0].plot(wlm, pca.flux_std + 1.5, "0.5")

ax[0].annotate(r"$\xi_\mu$", (2.02, 2.15))
ax[0].annotate(r"$\xi_\sigma$", (2.02, 1.6))

ax[0].set_ylabel(r"$\propto f_\lambda$ + offset")

for i, pcomp in enumerate(pca.pcomps):
    offset = (1.2 - i*0.25)
    ax[0].plot(wlm, pcomp + offset, "k")
    ax[0].annotate(r"${} = {:0.1f}$".format(w_list[i], ww[i]), (2.035, offset + 0.05), color="k", size=8)
    ax[0].annotate(xi_list[i], (2.02, offset + 0.05), color="k", size=9)

#ax[0].set_ylim(0.0, 2.15)
    
ax[1].plot(wlm, h_spec + 0.15, "b")
ax[1].plot(wlm, recon_spec, "g")
ax[1].annotate("synthetic", (2.02, 0.9), color="b", size=9)
ax[1].annotate("PCA reconstructed", (2.02, 0.8), color="g", size=9)
#ax[1].xaxis.set_major_formatter(FSF("%.1f"))
ax[1].set_xlabel(r"$\lambda$ [$\mu$m]")
ax[1].set_ylabel(r"$\propto f_\lambda$ + offset")

#Plots the residual spectrum
#ax[1].plot(wl, (h_spec - recon_spec) * 10)
fig.subplots_adjust(hspace=0.15, left=0.1, right=0.9, bottom=0.15, top=0.93)
ax[-1].set_xlim(wlm[0], wlm[-1])
saveall(fig, "../../plots/pca_reconstruct")
#plt.show()


    

In [ ]: