In [1]:

    

from Starfish.emulator import Emulator, PCAGrid
from Starfish.utils import saveall
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
from matplotlib.ticker import FormatStrFormatter as FSF
from matplotlib.ticker import MaxNLocator
from matplotlib.ticker import MultipleLocator

from Starfish.grid_tools import HDF5Interface


    

In [2]:

    

em = Emulator.open("../../libraries/PHOENIX_SPEX_M_PCA.hdf5") # All weights
pca = em.PCAGrid


    

In [3]:

    

temps = np.unique(pca.gparams[:,0])
#loggs = np.unique(pca.gparams[:,1])
#Zs = np.unique(pca.gparams[:,2])
#points = {"temp":temps, "logg":loggs, "Z":Zs}


    

In [4]:

    

# temps use just 2800 - 3200
logg = 5.0
Z = 0.0


    

In [5]:

    

int_temps = np.linspace(temps[0], temps[-1], num=40)
int_temps2 = np.linspace(3100, 3200, num=40)

weight_index = 4
weights = pca.w[weight_index]
pcomp = pca.pcomps[weight_index]
emw = em.WEs[weight_index] # The emulator for the first weight

nsamp = 50


    

In [13]:

    

mu, var = emw(np.array([3150, logg, Z]))
mu = mu[0]
var = var[0][0]


    

In [15]:

    

print(mu)


    

    




-2.19961725216

In [35]:

    

fig = plt.figure(figsize=(4, 3))
ax = fig.add_subplot(111)

# Load the grid-point weights
ww = []
for temp in temps:
    pars = np.array([temp, logg, Z])
    index = pca.get_index(pars)
    ww.append(weights[index])
    
fparams = []
for temp in int_temps:
    fparams.append([temp, logg, Z])
fparams = np.array(fparams)        

wgps = []
for k in range(nsamp):
    wgps.append(emw.draw_weights(fparams))
    
for k in range(nsamp):    
    ax.plot(int_temps, wgps[k] , "b", lw=0.1)
    
ax.plot(temps, ww, "bo")
    
ax.xaxis.set_major_formatter(FSF("%.0f"))        
ax.xaxis.set_major_locator(MultipleLocator(100))  

ax.set_xlabel(r"$T_\textrm{eff}$")
ax.set_ylabel(r"$w_5$")
        
#        #l,b,w,h
# rect = 0.6, 0.2, 0.25, 0.3

# axins = fig.add_axes(rect)

axins = inset_axes(ax, width=1.5, height=1.5, loc=4)

# Finer temperature spacing for the inset
fparams2 = []
for temp in int_temps2:
    fparams2.append([temp, logg, Z])
fparams2 = np.array(fparams2)       

for k in range(nsamp):    
    axins.plot(int_temps2, emw.draw_weights(fparams2) , "b", lw=0.1)

axins.axvline(3150, color="k", ymin=0.2, ymax=0.8)    
# sub region of the original image
x1, x2, y1, y2 = 3140, 3160, mu-0.01, mu+0.01
axins.set_xlim(x1, x2)
axins.set_ylim(y1, y2)
#ax.plot(3150, mu, "o")
axins.plot(3150, mu, "ko")

axins.xaxis.set_ticklabels([])
axins.yaxis.set_ticklabels([])

mark_inset(ax, axins, loc1=2, loc2=1, fc="none", ec="0.6")

fig.subplots_adjust(left=0.14, right=0.86, bottom=0.15)

saveall(fig, "../../plots/GP_left")

plt.show()


    

In [27]:

    

def gauss(x):
    return 1/np.sqrt(2. * np.pi * var) * np.exp(-0.5 * (x - mu)**2/var)


    

In [32]:

    

step = 0.1 * (y2 - y1)


    

In [33]:

    

xs = np.linspace(y1 + step, y2 - step, num=100)


    

In [34]:

    

fig = plt.figure(figsize=(3,3))
ax = fig.add_subplot(111)

ax.plot(xs, gauss(xs), "k")
ax.set_xlim(xs[0], xs[-1])

#ax.xaxis.set_major_formatter(FSF("%.0f"))        
ax.xaxis.set_major_locator(MultipleLocator(0.005))  
ax.set_xlabel(r"$w_5$")
ax.set_ylabel(r"$p(w_5 |\, \theta_\ast)$")

fig.subplots_adjust(left=0.18, right=0.82, bottom=0.15)

saveall(fig, "../../plots/GP_right")

plt.show()


    

In [ ]: