Plotting routines for paper



In [1]:

    
#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 [2]:

    
from StellarSpectra.model import Model
from StellarSpectra.spectrum import DataSpectrum
from StellarSpectra.grid_tools import TRES, HDF5Interface
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)
myOrderModel = myModel.OrderModels[0]
spec = myModel.get_data()

wl = spec.wls[0]
fl = spec.fls[0]

model_fl = myOrderModel.get_spectrum()

median_fl = np.median(fl)
fl_norm = fl / median_fl
model_fl_norm = model_fl / median_fl

residuals = fl_norm - model_fl_norm
cheb = myOrderModel.get_Cheb()









    



Determine Chunk Log: Wl is 8192
Creating OrderModel 0
params are  {'c3': -0.004788022228329032, 'c1': -0.017684987777693707, 'logc0': -0.030610219611502083, 'c2': -0.017530605981209466}

Poster plots of single spectra



In [29]:

    
figsize=(14, 4)
lw = 3

fig = plt.figure(figsize=figsize)
ax = fig.add_subplot(111)
ax.plot(wl, fl_norm, "b", lw=lw)

ax.set_xlim(5152, 5209)
#print(np.min(fl_norm), np.max(fl_norm))
ax.set_ylim(0.18, 1.09)
fig.subplots_adjust(left=0, right=1, bottom=0, top=1)

ax.axis('off')
fig.savefig("../plots/data.svg")
plt.show()



In [30]:

    
fig = plt.figure(figsize=figsize)
ax = fig.add_subplot(111)
ax.plot(wl, model_fl_norm, "r", lw=lw)
ax.set_ylim(0.18, 1.09)
ax.set_xlim(5152, 5209)

fig.subplots_adjust(left=0, right=1, bottom=0, top=1)

ax.axis('off')
fig.savefig("../plots/model.svg")
plt.show()



In [31]:

    
fig = plt.figure(figsize=figsize)
ax = fig.add_subplot(111)
ax.plot(wl, residuals, "k", lw=lw)
span = (1.09 - 0.18)/2
ax.set_ylim(-span, span)
ax.set_xlim(5152, 5209)

fig.subplots_adjust(left=0, right=1, bottom=0, top=1)

ax.axis('off')
fig.savefig("../plots/residuals.svg")
plt.show()

Large view residual plot



In [9]:

    
fig, ax = plt.subplots(nrows=2, figsize=(6,2.7), sharex=True)


l0, = ax[0].plot(wl, fl_norm, "b")
l1, = ax[0].plot(wl, model_fl_norm, "r")
ax[0].set_ylabel(r"$\propto f_\lambda$")
ax[0].yaxis.set_major_locator(MaxNLocator(nbins=7, prune='lower'))
ax[0].fill_betweenx(np.array([0.1, 1.2]), 5144, 5150, color="#FF8800", alpha=0.5)
ax[0].set_ylim(0.1, 1.2)


ax[1].fill_betweenx(np.array([-0.55, 0.2]), 5144, 5150, color="#FF8800", alpha=0.5)
l2, = ax[1].plot(wl, residuals, "k")
ax[1].xaxis.set_major_formatter(FSF("%.0f"))
ax[1].set_xlabel(r"$\lambda$ (\AA)")
ax[1].set_ylabel(r"residuals $\propto f_\lambda$")
ax[1].yaxis.set_major_locator(MaxNLocator(nbins=7, prune='upper'))
ax[1].set_xlim(5135, 5236)
ax[1].set_ylim(-0.55, 0.2)
ax[1].legend([l0, l1, l2], ["data", "model", "residuals"], loc="lower right", prop={'size':9})

fig.subplots_adjust(hspace=0, left=0.14, right=0.95, bottom=0.18, top=0.94)

fig.savefig("../plots/residual_23.png")
fig.savefig("../plots/residual_23.pdf")
fig.savefig("../plots/residual_23.svg")
plt.show()

Proposal Plot of Data + Model, no residuals

Chebyshev, Data + Model, and Residuals



In [29]:

    
nullfmt = matplotlib.ticker.NullFormatter()
import matplotlib as mpl
mpl.rc('lines', linewidth=1)

fig = plt.figure(figsize=(7,5))
#[left, bottom, width, height]
left = 0.12
width = 0.85

plt_base = 0.13
label_ax = fig.add_axes([left - 0.1, 0.05, .95, .95],alpha=0.0,frame_on=False,xticks=[],yticks=[])

#unit of height for each plot
height = 0.15

ax_resid = fig.add_axes([left,plt_base,width, 1.3 * height])
ax_resid.plot(wl, residuals, "k", label="residuals")
ax_resid.xaxis.set_major_formatter(FSF("%.0f"))
ax_resid.set_xlabel(r"$\lambda$  (\AA)")
#ax_resid.legend(loc='lower left')
ax_resid.text(0.02, 0.1, "residuals", transform=ax_resid.transAxes)
ax_resid.set_ylabel("$\propto f_\lambda$")
ax_resid.yaxis.set_major_locator(MaxNLocator(nbins=5, prune='upper'))

ax_data = fig.add_axes([left,plt_base + 1.3 * height,width, 3 * height])
ax_data.xaxis.set_major_formatter(nullfmt)
ax_data.plot(wl, fl_norm, "b", label="data")
ax_data.plot(wl, model_fl_norm, "r:", label="model")
lg = ax_data.legend(loc="lower left")
lg.draw_frame(False)
ax_data.set_ylabel("$\propto f_\lambda$")
ax_data.yaxis.set_major_locator(MaxNLocator(prune='both'))

ax_cheb = fig.add_axes([left,plt_base + 4.3 * height,width, height])
ax_cheb.xaxis.set_major_formatter(nullfmt)
ax_cheb.plot(wl, cheb, "k", label="Chebyshev")
ax_cheb.set_ylim(.99, 1.02)
#ax_cheb.legend(loc="lower left")
ax_cheb.text(0.02, 0.1, "Chebyshev", transform=ax_cheb.transAxes)
ax_cheb.yaxis.set_major_locator(MaxNLocator(nbins=5, prune='lower'))

for ax in [ax_resid, ax_data, ax_cheb]:
    ax.set_xlim(5155,5200)

fig.savefig("../plots/model_data.eps")
plt.show()

Examples of correlated noise and bad regions



In [13]:

    
import matplotlib as mpl
mpl.rc('lines', linewidth=1)

fig, ax = plt.subplots(nrows=2, ncols=4, figsize=(6.5,4.))
spec2 = DataSpectrum.open("WASP14/WASP-14_2009-06-15_04h13m57s_cb.spec.flux", orders=[22,23])
fl23, fl24 = spec2.fls

median_fl23 = np.median(fl23)
median_fl24 = np.median(fl24)

fl23n = fl23/median_fl23
fl24n = fl24/median_fl24

#Model
wl23 = np.load("../residuals/wlsz23.npy")[0]
f23 = np.load("../residuals/fl23.npy")[0]/median_fl23
wl24 = np.load("../residuals/wlsz24.npy")[0]
f24 = np.load("../residuals/fl24.npy")[0]/median_fl24


r23 = (fl23n - f23)
r24 = (fl24n - f24)

ax[0,0].plot(wl23, fl23n, "b", label="data", drawstyle="steps-mid")
ax[0,0].plot(wl23, f23, "r:", label="model", drawstyle="steps-mid")

ax[1,0].plot(wl23, r23, "k", drawstyle="steps-mid")
#ax[0,0].set_xlim(5136.4, 5139)
#ax[1,0].set_xlim(5136.4, 5139)
ax[1,0].set_ylim(-.1,.1)

ax[0,0].set_xlim(5144, 5150)
ax[1,0].set_xlim(5144, 5150)

ax[0,1].plot(wl23, fl23n, "b", drawstyle="steps-mid")
ax[0,1].plot(wl23, f23, "r:", drawstyle="steps-mid")
ax[1,1].plot(wl23, r23, "k", drawstyle="steps-mid")
ax[0,1].set_xlim(5188, 5189.5)
ax[1,1].set_xlim(5188, 5189.5)
#ax[1,1].set_ylim(-4,4)

ax[0,2].plot(wl24, fl24n, "b", label='data', drawstyle="steps-mid")
ax[0,2].plot(wl24, f24, "r:", label='model', drawstyle="steps-mid")
ax[0,2].legend(loc="lower center", prop={'size':10})

ax[1,2].plot(wl24, r24, "k", drawstyle="steps-mid")
ax[0,2].set_xlim(5258, 5260)
ax[1,2].set_xlim(5258, 5260)

ax[0,3].plot(wl24, fl24n, "b", drawstyle="steps-mid")
ax[0,3].plot(wl24, f24, "r:", drawstyle="steps-mid")
ax[1,3].plot(wl24, r24, "k", drawstyle="steps-mid")
ax[0,3].set_xlim(5260, 5269)
ax[1,3].set_xlim(5260, 5269)
#ax[0,3].set_xlim(5260, 5271)
#ax[1,3].set_xlim(5260, 5271)
ax[1,3].set_ylim(-.18, .18)


for j in range(4):
    labels = ax[1,j].get_xticklabels()
    for label in labels:
        label.set_rotation(60)

ax[0,0].set_ylabel(r"$\propto f_\lambda$")
ax[1,0].set_ylabel(r"residuals $\propto f_\lambda$")

for i in range(2):
    for j in range(4):
        ax[i,j].xaxis.set_major_formatter(FSF("%.0f"))
        ax[i,j].xaxis.set_major_locator(MultipleLocator(1.))
        ax[i,j].tick_params(axis='both', which='major', labelsize=10)

for i in [0,1]:
    for j in [1,2]:
        ax[i,j].xaxis.set_major_formatter(FSF("%.1f"))
        ax[i,j].xaxis.set_major_locator(MultipleLocator(0.5))

for i in [0,1]:
    ax[i,3].xaxis.set_major_formatter(FSF("%.0f"))
    ax[i,3].xaxis.set_major_locator(MultipleLocator(2))

class_label = ["0", "I", "II", "III"]
for j in range(4):
    ax[0,j].set_title(class_label[j])
    ax[0,j].xaxis.set_ticklabels([])
    ax[0,j].set_ylim(0.4, 1.15)
    if j != 0:
        ax[0,j].yaxis.set_ticklabels([])

fig.subplots_adjust(left=0.13, right=0.99, top=0.94, bottom=0.18, hspace=0.1, wspace=0.3)
fig.text(0.48, 0.02, r"$\lambda$ (\AA)")
fig.savefig("../plots/badlines.pdf")
plt.show()

Autocorrelation of Class 0 region

Create two separate figures, such that they can be shown side-by-side using \includegraphics{}



In [10]:

    
fig, ax = plt.subplots(nrows=2, sharex=True, figsize=(3.5, 3.5))
spec2 = DataSpectrum.open("../data/WASP14/WASP14-2009-06-14.hdf5", orders=[22])
fl23 = spec2.fls[0]

median_fl23 = np.median(fl23)

fl23n = fl23/median_fl23

#Model
wl23 = np.load("../residuals/wlsz23.npy")[0]
f23 = np.load("../residuals/fl23.npy")[0]/median_fl23

ind = (wl23 > 5144) & (wl23 < 5150)
wl23 = wl23[ind]
f23 = f23[ind]
fl23n = fl23n[ind]

r23 = (fl23n - f23)

ax[0].plot(wl23, fl23n, "b", label="data", drawstyle="steps-mid")
ax[0].plot(wl23, f23, "r", label="model", drawstyle="steps-mid")
#ax[0].legend(loc="lower left")
ax[0].set_ylabel(r"$\propto f_\lambda$")


ax[1].plot(wl23, r23, "k", label="residuals", drawstyle="steps-mid")
#ax[1].legend(loc="lower right")
ax[1].set_ylabel(r"residuals $\propto f_\lambda$")
ax[1].set_xlabel(r"$\lambda$ [\AA]")

ax[1].xaxis.set_major_formatter(FSF("%.0f"))

fig.subplots_adjust(right=0.92, bottom=0.13, left=0.22)
fig.savefig("../plots/class0_residuals.png")
fig.savefig("../plots/class0_residuals.svg")
fig.savefig("../plots/class0_residuals.pdf")
plt.show()

$$\hat{R}(k) = \frac{1}{(n - k) \sigma^2} \sum_{t = 1}^{n - k} (X_t - \mu) (X_{t+k} - \mu) $$



In [63]:

    
def estimated_autocorrelation(x):
    '''From http://stackoverflow.com/questions/14297012/estimate-autocorrelation-using-python'''
    n = len(x)
    variance = np.var(x)
    x = x - np.mean(x)
    r = np.correlate(x, x, mode = 'full')[-n:] #takes the positive half of the correlation
    #assert N.allclose(r, N.array([(x[:n-k]*x[-(n-k):]).sum() for k in range(n)]))
    result = r/(variance*(np.arange(n, 0, -1))) #divides by the number pixels used in the estimation
    return result



In [67]:

    
#Now estimate the autocorrelation of the residuals 
fig = plt.figure(figsize=(2.5,2.5))
ax = fig.add_subplot(111)
ax.plot(estimated_autocorrelation(r23), "m", drawstyle="steps-mid")
ax.set_xlabel(r"pixel offset $k$")
ax.set_ylabel(r"$\hat{\rho}(k)$")
ax.set_xlim(-1,60)
fig.subplots_adjust(left=0.27, bottom=0.17, right=0.95)
fig.savefig("../plots/class0_autocorrelation.png")
fig.savefig("../plots/class0_autocorrelation.pdf")
fig.savefig("../plots/class0_autocorrelation.svg")
plt.show()

Global covariance matrix



In [32]:

    
def k_3_2(r, l):
    return (1 + np.sqrt(3) * r/l) * np.exp(-np.sqrt(3) * r/l)

@np.vectorize
def hann(r, r0):
    if np.abs(r) < r0:
        return 0.5 + 0.5 * np.cos(np.pi * r/r0)
    else:
        return 0

def k_3_2_hann(r, l):
    return k_3_2(r, l) * hann(r, 6 * l)



In [33]:

    
def gauss_func(x0i, x1i, x0v=None, x1v=None, amp=None, mu=None, sigma=None):
    x0 = x0v[x0i]
    x1 = x1v[x1i]
    return amp**2/(2 * np.pi * sigma**2) * np.exp(-((x0 - mu)**2 + (x1 - mu)**2)/(2 * sigma**2))

def k_3_2_func(x0i, x1i, x0v=None, x1v=None, amp=None, l=None):
    x0 = x0v[x0i]
    x1 = x1v[x1i]
    r = np.abs(x1 - x0)
    return amp**2 * k_3_2(r, l)

def k_3_2_hann_func(x0i, x1i, x0v=None, x1v=None, amp=None, l=None):
    x0 = x0v[x0i]
    x1 = x1v[x1i]
    r = np.abs(x1 - x0)
    return amp**2 * k_3_2_hann(r, l)

Matrix Thumbnail Plot for poster



In [37]:

    
spec2 = DataSpectrum.open("WASP14/WASP-14_2009-06-15_04h13m57s_cb.spec.flux", orders=[22])
fl23 = spec2.fls[0]
sigma23 = spec2.sigmas[0]

median_fl23 = np.median(fl23)

fl23n = fl23/median_fl23
sigma23n = sigma23/median_fl23

#Model
wl23 = np.load("../residuals/wlsz23.npy")[0]
f23 = np.load("../residuals/fl23.npy")[0]/median_fl23

#ind = (wl23 > 5144) & (wl23 < 5150)
ind = (wl23 > 5144) & (wl23 < 5145)
wl23 = wl23[ind]
f23 = f23[ind]
fl23n = fl23n[ind]
sigma23n = sigma23n[ind]

r23 = (fl23n - f23)

N = len(wl23)



In [38]:

    
mat_3_2 = np.fromfunction(k_3_2_hann_func, (N,N), x0v=wl23, x1v=wl23, amp=0.028, l=0.14, dtype=np.int) #matrix from Matern kernel
mat = mat_3_2 + sigma23n**2 * np.eye(N) #add in the per-pixel noise



In [55]:

    
fig = plt.figure(figsize=(14, 14))
ax = fig.add_subplot(111)
mi, ma = wl23[0], wl23[-1]
img = ax.imshow(mat, origin="upper", interpolation="none", extent=[mi, ma, ma, mi], cmap=plt.cm.OrRd)

cax = fig.add_axes([0.93, 0.2, 0.06, 0.6])
cb = fig.colorbar(img, cax=cax)
cb.set_ticks([])

fig.subplots_adjust(left=0.02, bottom=0.02, top=0.98, right=0.88)
ax.axis('off')
fig.savefig("../plots/matrix_thumbnail.pdf")
#plt.show()

Matrix full plot for poster



In [57]:

    
spec2 = DataSpectrum.open("WASP14/WASP-14_2009-06-15_04h13m57s_cb.spec.flux", orders=[22])
fl23 = spec2.fls[0]
sigma23 = spec2.sigmas[0]

median_fl23 = np.median(fl23)

fl23n = fl23/median_fl23
sigma23n = sigma23/median_fl23

#Model
wl23 = np.load("../residuals/wlsz23.npy")[0]
f23 = np.load("../residuals/fl23.npy")[0]/median_fl23


ind = (wl23 > 5152) & (wl23 < 5209)
wl23 = wl23[ind]
f23 = f23[ind]
fl23n = fl23n[ind]
sigma23n = sigma23n[ind]

r23 = (fl23n - f23)

N = len(wl23)



In [58]:

    
mat_3_2 = np.fromfunction(k_3_2_hann_func, (N,N), x0v=wl23, x1v=wl23, amp=0.028, l=0.14, dtype=np.int) #matrix from Matern kernel
mat = mat_3_2 + sigma23n**2 * np.eye(N) #add in the per-pixel noise



In [59]:

    
fig = plt.figure(figsize=(14, 14))
ax = fig.add_subplot(111)
mi, ma = wl23[0], wl23[-1]
img = ax.imshow(mat, origin="upper", interpolation="none", extent=[mi, ma, ma, mi], cmap=plt.cm.OrRd)

cax = fig.add_axes([0.93, 0.2, 0.06, 0.6])
cb = fig.colorbar(img, cax=cax)
cb.set_ticks([])

fig.subplots_adjust(left=0.02, bottom=0.02, top=0.98, right=0.88)
ax.axis('off')
fig.savefig("../plots/matrix_full.pdf")
#plt.show()

Plot for paper

Using the same wavelength range from the Class 0 plot, imshow the Matern matrix



In [34]:

    
spec2 = DataSpectrum.open("WASP14/WASP-14_2009-06-15_04h13m57s_cb.spec.flux", orders=[22])
fl23 = spec2.fls[0]
sigma23 = spec2.sigmas[0]

median_fl23 = np.median(fl23)

fl23n = fl23/median_fl23
sigma23n = sigma23/median_fl23

#Model
wl23 = np.load("../residuals/wlsz23.npy")[0]
f23 = np.load("../residuals/fl23.npy")[0]/median_fl23

#ind = (wl23 > 5144) & (wl23 < 5150)
ind = (wl23 > 5144) & (wl23 < 5145)
wl23 = wl23[ind]
f23 = f23[ind]
fl23n = fl23n[ind]
sigma23n = sigma23n[ind]

r23 = (fl23n - f23)

N = len(wl23)



In [6]:

    
mat_3_2 = np.fromfunction(k_3_2_hann_func, (N,N), x0v=wl23, x1v=wl23, amp=0.028, l=0.14, dtype=np.int) #matrix from Matern kernel
mat = mat_3_2 + sigma23n**2 * np.eye(N) #add in the per-pixel noise



In [18]:

    
print(np.max(mat))









    



0.00108466625491



In [19]:

    
fig = plt.figure(figsize=(3.5,3.5))
ax = fig.add_subplot(111)
mi, ma = wl23[0], wl23[-1]
img = ax.imshow(mat, origin="upper", interpolation="none", extent=[mi, ma, ma, mi], cmap=plt.cm.OrRd)
ax.xaxis.set_major_formatter(FSF("%.1f"))

labels = ax.get_xticklabels()
for label in labels:
    label.set_rotation(40)

ax.yaxis.set_major_formatter(FSF("%.1f"))

ax.set_xlabel(r"$\lambda$ [\AA]")
ax.set_ylabel(r"$\lambda$ [\AA]")

cax = fig.add_axes([0.85, 0.2, 0.03, 0.65])
cb = fig.colorbar(img, cax=cax)

ticks = np.linspace(0, 1e-3, num=6)
cb.set_ticks(ticks)

fig.suptitle(r"$\sigma_{ij}$",x=0.76,y=0.89)

fig.subplots_adjust(left=0.24, bottom=0.13, top=0.97, right=0.8)
fig.savefig("../plots/matern_matrix.pdf")
fig.savefig("../plots/matern_matrix.png")
plt.show()

Draws from the Global covariance matrix



In [20]:

    
spec2 = DataSpectrum.open("WASP14/WASP-14_2009-06-15_04h13m57s_cb.spec.flux", orders=[22])
fl23 = spec2.fls[0]
sigma23 = spec2.sigmas[0]

median_fl23 = np.median(fl23)

fl23n = fl23/median_fl23
sigma23n = sigma23/median_fl23

#Model
wl23 = np.load("../residuals/wlsz23.npy")[0]
f23 = np.load("../residuals/fl23.npy")[0]/median_fl23

ind = (wl23 > 5144) & (wl23 < 5150)
wl23 = wl23[ind]
f23 = f23[ind]
fl23n = fl23n[ind]
sigma23n = sigma23n[ind]

r23 = (fl23n - f23)

N = len(wl23)



In [41]:

    
mat_3_2 = np.fromfunction(k_3_2_hann_func, (N,N), x0v=wl23, x1v=wl23, amp=0.028, l=0.14, dtype=np.int) #matrix from Matern kernel
mat = mat_3_2 + sigma23n**2 * np.eye(N) #add in the per-pixel noise



In [52]:

    
fl_fake = np.random.multivariate_normal(np.zeros((N,)), mat)
fl_fake2 = np.random.multivariate_normal(np.zeros((N,)), mat)
fl_fake3 = np.random.multivariate_normal(np.zeros((N,)), mat)



In [53]:

    
fig, ax = plt.subplots(nrows=4, sharex=True, figsize=(3.2, 3.2))

left = 5144.1
ax[0].plot(wl23, r23, "k", label="residuals", drawstyle="steps-mid")
ax[0].text(left, 0.04, "residuals")

bot = -0.09
ax[1].plot(wl23, fl_fake1, "b", drawstyle="steps-mid", label="random draw")
ax[1].text(left, bot, "random draw")

ax[2].plot(wl23, fl_fake2, "b", drawstyle="steps-mid")
ax[2].text(left, bot, "random draw")

ax[3].plot(wl23, fl_fake3, "b", drawstyle="steps-mid")
ax[3].text(left, bot, "random draw")

#tl = [text.get_text() for text in ax[0].yaxis.get_ticklabels()]
#print(tl)

for j in range(3):
    ax[j + 1].yaxis.set_ticklabels([])
    
for j in range(4):
    ax[j].yaxis.set_major_locator(MaxNLocator(nbins=5))
    ax[j].set_ylim(-0.1,0.1)
    

ax[2].set_ylabel(r"$\propto f_\lambda$")
ax[-1].set_xlabel(r"$\lambda$ [\AA]")
ax[-1].xaxis.set_major_formatter(FSF("%.0f"))
ax[-1].xaxis.set_major_locator(MultipleLocator(1.))
ax[-1].set_xlim(5144, 5149.5)

fig.subplots_adjust(hspace=0., right=0.98, top=0.95, bottom=0.15, left=0.15)
fig.savefig("../plots/matern_draw.pdf")
plt.show()

Global covariance movie



In [5]:

    
spec2 = DataSpectrum.open("WASP14/WASP-14_2009-06-15_04h13m57s_cb.spec.flux", orders=[22])
fl23 = spec2.fls[0]
sigma23 = spec2.sigmas[0]

median_fl23 = np.median(fl23)

fl23n = fl23/median_fl23
sigma23n = sigma23/median_fl23

#Model
wl23 = np.load("../residuals/wlsz23.npy")[0]
f23 = np.load("../residuals/fl23.npy")[0]/median_fl23

ind = (wl23 > 5144) & (wl23 < 5150)
wl23 = wl23[ind]
f23 = f23[ind]
fl23n = fl23n[ind]
sigma23n = sigma23n[ind]

r23 = (fl23n - f23)

N = len(wl23)



In [6]:

    
mat_3_2 = np.fromfunction(k_3_2_hann_func, (N,N), x0v=wl23, x1v=wl23, amp=0.028, l=0.14, dtype=np.int) #matrix from Matern kernel
mat = mat_3_2 + sigma23n**2 * np.eye(N) #add in the per-pixel noise



In [30]:

    
# First set up the figure, the axis, and the plot element we want to animate
fig, ax = plt.subplots(nrows=2, figsize=(4,4), sharex=True, sharey=True)

ax[0].plot(wl23, r23, "k", label="residuals", drawstyle="steps-mid")
ax[0].set_title("data residuals")


#Function to generate interpolated spectra
container = np.empty((2, len(wl23)))
def interp_spec(spec0, spec1, percentage):
    '''
    Return the interpolated spectrum between two spectra, weighted by
    the percentage (out of 1) of elapsed time between spectra.
    '''
    container[0, :] = spec0 * (1 - percentage)
    container[1, :] = spec1 * percentage
    return np.sum(container, axis=0)

ndraws = 10 #How many random spectra to draw
trans = np.concatenate((np.zeros((50,)), np.linspace(0, 1, num=5)))
print(trans)
ntrans = len(trans) 

#Create an array to hold all of the spectra
#starts with 1 + (ndraws - 1) * ntrans 
shape = (1 + ndraws * ntrans, len(wl23))
all_spectra = np.empty(shape)
#Fill it with fake draws and transitions between them

fake0 = np.random.multivariate_normal(np.zeros((N,)), mat)
all_spectra[0, :] = fake0
counter = 1
for draw in range(ndraws):
    fake1 = np.random.multivariate_normal(np.zeros((N,)), mat)
    for tran in trans:
        all_spectra[counter, :] = interp_spec(fake0, fake1, tran)
        counter += 1
    fake0 = fake1
        
#fl_fake = np.random.multivariate_normal(np.zeros((N,)), mat)
line, = ax[1].plot(wl23, fl_fakes[0], "b", drawstyle="steps-mid")


ax[1].set_title("fake residuals")
ax[-1].set_xlabel(r"$\lambda$ [\AA]")
ax[-1].xaxis.set_major_formatter(FSF("%.0f"))
ax[-1].xaxis.set_major_locator(MultipleLocator(1.))
ax[-1].set_ylim(-0.1, 0.1)
ax[-1].set_xlim(5144.5, 5149.5)
fig.subplots_adjust(top=0.93, bottom=0.12, right=0.97)

for i,spec in enumerate(all_spectra[1:]):
    line.set_ydata(spec)
    
    fig.canvas.draw_idle()
    fig.savefig("../plots/animation/global/{:0>3d}.png".format(i))
    #plt.close(fig)









    



[ 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
  0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
  0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
  0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
  0.25  0.5   0.75  1.  ]
shape is (551, 128)
draw # 0
draw # 1
draw # 2
draw # 3
draw # 4
draw # 5
draw # 6
draw # 7
draw # 8
draw # 9



In [31]:

    
#Generating movies using ffmpeg
!ffmpeg -y -r 30 -i ../plots/animation/global/%03d.png -vcodec libx264 -r 30 -pix_fmt yuv420p ../plots/animation/global/global_animation.mp4









    



ffmpeg version 2.2.1 Copyright (c) 2000-2014 the FFmpeg developers
  built on Apr 10 2014 16:26:13 with gcc 4.8.2 (GCC) 20140206 (prerelease)
  configuration: --prefix=/usr --disable-debug --disable-static --enable-avresample --enable-dxva2 --enable-fontconfig --enable-gnutls --enable-gpl --enable-libass --enable-libbluray --enable-libfreetype --enable-libgsm --enable-libmodplug --enable-libmp3lame --enable-libopencore_amrnb --enable-libopencore_amrwb --enable-libopenjpeg --enable-libopus --enable-libpulse --enable-librtmp --enable-libschroedinger --enable-libspeex --enable-libtheora --enable-libv4l2 --enable-libvorbis --enable-libvpx --enable-libx264 --enable-libx265 --enable-libxvid --enable-pic --enable-postproc --enable-runtime-cpudetect --enable-shared --enable-swresample --enable-vdpau --enable-version3 --enable-x11grab
  libavutil      52. 66.100 / 52. 66.100
  libavcodec     55. 52.102 / 55. 52.102
  libavformat    55. 33.100 / 55. 33.100
  libavdevice    55. 10.100 / 55. 10.100
  libavfilter     4.  2.100 /  4.  2.100
  libavresample   1.  2.  0 /  1.  2.  0
  libswscale      2.  5.102 /  2.  5.102
  libswresample   0. 18.100 /  0. 18.100
  libpostproc    52.  3.100 / 52.  3.100
Input #0, image2, from '../plots/animation/%03d.png':
  Duration: 00:00:22.00, start: 0.000000, bitrate: N/A
    Stream #0:0: Video: png, rgba, 800x800 [SAR 7874:7874 DAR 1:1], 25 fps, 25 tbr, 25 tbn, 25 tbc
[libx264 @ 0x11db500] using SAR=1/1
[libx264 @ 0x11db500] using cpu capabilities: MMX2 SSE2Fast SSSE3 SSE4.2
[libx264 @ 0x11db500] profile High, level 3.1
[libx264 @ 0x11db500] 264 - core 142 r2397 aff928d - H.264/MPEG-4 AVC codec - Copyleft 2003-2014 - http://www.videolan.org/x264.html - options: cabac=1 ref=3 deblock=1:0:0 analyse=0x3:0x113 me=hex subme=7 psy=1 psy_rd=1.00:0.00 mixed_ref=1 me_range=16 chroma_me=1 trellis=1 8x8dct=1 cqm=0 deadzone=21,11 fast_pskip=1 chroma_qp_offset=-2 threads=6 lookahead_threads=1 sliced_threads=0 nr=0 decimate=1 interlaced=0 bluray_compat=0 constrained_intra=0 bframes=3 b_pyramid=2 b_adapt=1 b_bias=0 direct=1 weightb=1 open_gop=0 weightp=2 keyint=250 keyint_min=25 scenecut=40 intra_refresh=0 rc_lookahead=40 rc=crf mbtree=1 crf=23.0 qcomp=0.60 qpmin=0 qpmax=69 qpstep=4 ip_ratio=1.40 aq=1:1.00
Output #0, mp4, to '../plots/animation/out.mp4':
  Metadata:
    encoder         : Lavf55.33.100
    Stream #0:0: Video: h264 (libx264) ([33][0][0][0] / 0x0021), yuv420p, 800x800 [SAR 1:1 DAR 1:1], q=-1--1, 15360 tbn, 30 tbc
Stream mapping:
  Stream #0:0 -> #0:0 (png -> libx264)
Press [q] to stop, [?] for help

video:244kB audio:0kB subtitle:0 data:0 global headers:0kB muxing overhead 2.946554%
[libx264 @ 0x11db500] frame I:3     Avg QP:19.37  size: 28076
[libx264 @ 0x11db500] frame P:149   Avg QP:15.08  size:   776
[libx264 @ 0x11db500] frame B:398   Avg QP:14.62  size:   123
[libx264 @ 0x11db500] consecutive B-frames:  3.3%  0.4%  1.1% 95.3%
[libx264 @ 0x11db500] mb I  I16..4: 31.5% 42.7% 25.8%
[libx264 @ 0x11db500] mb P  I16..4:  0.4%  0.3%  0.4%  P16..4:  0.7%  0.3%  0.1%  0.0%  0.0%    skip:97.8%
[libx264 @ 0x11db500] mb B  I16..4:  0.1%  0.0%  0.0%  B16..8:  1.5%  0.1%  0.0%  direct: 0.0%  skip:98.4%  L0:59.2% L1:40.6% BI: 0.1%
[libx264 @ 0x11db500] 8x8 transform intra:35.8% inter:12.6%
[libx264 @ 0x11db500] coded y,uvDC,uvAC intra: 14.1% 19.5% 18.3% inter: 0.1% 0.2% 0.2%
[libx264 @ 0x11db500] i16 v,h,dc,p: 77% 21%  2%  0%
[libx264 @ 0x11db500] i8 v,h,dc,ddl,ddr,vr,hd,vl,hu: 54% 13% 32%  0%  0%  0%  0%  0%  0%
[libx264 @ 0x11db500] i4 v,h,dc,ddl,ddr,vr,hd,vl,hu: 59% 17% 15%  1%  2%  2%  2%  1%  1%
[libx264 @ 0x11db500] i8c dc,h,v,p: 74%  4% 22%  0%
[libx264 @ 0x11db500] Weighted P-Frames: Y:0.0% UV:0.0%
[libx264 @ 0x11db500] ref P L0: 76.1%  3.3% 15.7%  5.0%
[libx264 @ 0x11db500] ref B L0: 59.1% 39.3%  1.6%
[libx264 @ 0x11db500] ref B L1: 98.1%  1.9%
[libx264 @ 0x11db500] kb/s:108.52

Gaussian line kernel



In [3]:

    
def k_3_2(r, l):
    return (1 + np.sqrt(3) * r/l) * np.exp(-np.sqrt(3) * r/l)

@np.vectorize
def hann(r, r0):
    if np.abs(r) < r0:
        return 0.5 + 0.5 * np.cos(np.pi * r/r0)
    else:
        return 0

def k_3_2_hann(r, l):
    return k_3_2(r, l) * hann(r, 6 * l)



In [4]:

    
def gauss_func(x0i, x1i, x0v=None, x1v=None, amp=None, mu=None, sigma=None):
    x0 = x0v[x0i]
    x1 = x1v[x1i]
    return amp**2/(2 * np.pi * sigma**2) * np.exp(-((x0 - mu)**2 + (x1 - mu)**2)/(2 * sigma**2))

def k_3_2_func(x0i, x1i, x0v=None, x1v=None, amp=None, l=None):
    x0 = x0v[x0i]
    x1 = x1v[x1i]
    r = np.abs(x1 - x0)
    return amp**2 * k_3_2(r, l)

def k_3_2_hann_func(x0i, x1i, x0v=None, x1v=None, amp=None, l=None):
    x0 = x0v[x0i]
    x1 = x1v[x1i]
    r = np.abs(x1 - x0)
    return amp**2 * k_3_2_hann(r, l)

def gauss_hann_func(x0i, x1i, x0v=None, x1v=None, amp=None, mu=None, sigma=None):
    x0 = x0v[x0i]
    x1 = x1v[x1i]
    r = np.abs(x1 - x0)
    return hann(r, 4 * sigma) * amp**2/(2 * np.pi * sigma**2) * np.exp(-((x0 - mu)**2 + (x1 - mu)**2)/(2 * sigma**2))



In [5]:

    
spec2 = DataSpectrum.open("WASP14/WASP-14_2009-06-15_04h13m57s_cb.spec.flux", orders=[22])
fl23 = spec2.fls[0]
sigma23 = spec2.sigmas[0]

median_fl23 = np.median(fl23)

fl23n = fl23/median_fl23
sigma23n = sigma23/median_fl23

#Model
wl23 = np.load("../residuals/wlsz23.npy")[0]
f23 = np.load("../residuals/fl23.npy")[0]/median_fl23

ind = (wl23 > 5188.1) & (wl23 < 5190.2)
wl23 = wl23[ind]
f23 = f23[ind]
fl23n = fl23n[ind]
sigma23n = sigma23n[ind]

r23 = (fl23n - f23)

N = len(wl23)



In [6]:

    
mat_3_2 = np.fromfunction(k_3_2_hann_func, (N,N), x0v=wl23, x1v=wl23, amp=0.028, l=0.14, dtype=np.int) #matrix from Matern kernel
mat_gauss = np.fromfunction(gauss_hann_func, (N,N), x0v=wl23, x1v=wl23, amp=0.09, mu=5188.84, sigma=0.1, dtype=np.int) #matrix from Matern kernel
mat = mat_3_2 + mat_gauss + sigma23n**2 * np.eye(N) #add in the per-pixel noise



In [18]:

    
from matplotlib.colors import LogNorm
fig = plt.figure(figsize=(3.5,3.5))
ax = fig.add_subplot(111)
mi, ma = wl23[0], wl23[-1]
img = ax.imshow(mat + 1e-4, origin="upper", interpolation="none", extent=[mi, ma, ma, mi], cmap=plt.cm.OrRd, norm=LogNorm(vmin=0.0001, clip=True))
ax.xaxis.set_major_formatter(FSF("%.1f"))

labels = ax.get_xticklabels()
for label in labels:
    label.set_rotation(40)

ax.yaxis.set_major_formatter(FSF("%.1f"))


ax.set_xlabel(r"$\lambda$ [\AA]")
ax.set_ylabel(r"$\lambda$ [\AA]")

cax = fig.add_axes([0.85, 0.2, 0.03, 0.65])
cb = fig.colorbar(img, cax=cax)

#Set two colorbars in an image? 
#ticks = np.linspace(0, 1e-3, num=6)
#cb.set_ticks(ticks)

fig.suptitle(r"$\sigma_{ij}$",x=0.76,y=0.9)

fig.subplots_adjust(left=0.24, bottom=0.13, top=0.97, right=0.8)
fig.savefig("../plots/gauss_matrix.pdf")
fig.savefig("../plots/gauss_matrix.png")
plt.show()

Draws from the Gaussian line kernel



In [17]:

    
fig = plt.figure()
ax = fig.add_subplot(111)

for i in range(20):
    fl_fake = np.random.multivariate_normal(np.zeros((N,)), mat)
    ax.plot(wl23, fl_fake, "k", lw=0.2)


ax.plot(wl23, r23)
plt.show()



In [23]:

    
gauss = 0.09 / (np.sqrt(2 * np.pi) * 0.1) * np.exp(-0.5 * (wl23 - 5188.8)**2/(0.01))
plt.plot(wl23, gauss)
plt.show()



In [37]:

    
fl_fake1 = np.random.multivariate_normal(np.zeros((N,)), mat)



In [38]:

    
fl_fake2 = np.random.multivariate_normal(np.zeros((N,)), mat)



In [28]:

    
fl_fake3 = np.random.multivariate_normal(np.zeros((N,)), mat)



In [44]:

    
fig, ax = plt.subplots(nrows=4, sharex=True, figsize=(3.2, 3.2))

right = 5189.4
top = -0.5
ax[0].plot(wl23, r23, "k", label="residuals", drawstyle="steps-mid")
ax[0].text(right, top, "residuals")


ax[1].plot(wl23, fl_fake1, "b", drawstyle="steps-mid", label="random draw")
ax[1].text(right, top, "random draw")

ax[2].plot(wl23, fl_fake2, "b", drawstyle="steps-mid")
ax[2].text(right, top, "random draw")

ax[3].plot(wl23, fl_fake3, "b", drawstyle="steps-mid")
ax[3].text(right, top, "random draw")

for j in range(3):
    ax[j + 1].yaxis.set_ticklabels([])
    
for j in range(4):
    ax[j].yaxis.set_major_locator(MaxNLocator(nbins=5))
    ax[j].set_ylim(-0.6,0.3)
    

ax[2].set_ylabel(r"$\propto f_\lambda$")
ax[-1].set_xlabel(r"$\lambda$ [\AA]")
ax[-1].xaxis.set_major_formatter(FSF("%.1f"))
ax[-1].xaxis.set_major_locator(MultipleLocator(0.5))
ax[-1].set_xlim(5188.1, 5190.2)

fig.subplots_adjust(hspace=0., right=0.98, top=0.95, bottom=0.15, left=0.15)
fig.savefig("../plots/gauss_draw.pdf")
fig.savefig("../plots/gauss_draw.png")
plt.show()

Gaussian line movie



In [52]:

    
# First set up the figure, the axis, and the plot element we want to animate
fig, ax = plt.subplots(nrows=2, figsize=(4,4), sharex=True, sharey=True)

ax[0].plot(wl23, r23, "k", label="residuals", drawstyle="steps-mid")
ax[0].set_title("data residuals")


#Function to generate interpolated spectra
container = np.empty((2, len(wl23)))
def interp_spec(spec0, spec1, percentage):
    '''
    Return the interpolated spectrum between two spectra, weighted by
    the percentage (out of 1) of elapsed time between spectra.
    '''
    container[0, :] = spec0 * (1 - percentage)
    container[1, :] = spec1 * percentage
    return np.sum(container, axis=0)

ndraws = 15 #How many random spectra to draw
trans = np.concatenate((np.zeros((15,)), np.linspace(0, 1, num=5)))
print(trans)
ntrans = len(trans) 

#Create an array to hold all of the spectra
#starts with 1 + (ndraws - 1) * ntrans 
shape = (1 + ndraws * ntrans, len(wl23))
print("shape is", shape)
all_spectra = np.empty(shape)
#Fill it with fake draws and transitions between them

fake0 = np.random.multivariate_normal(np.zeros((N,)), mat)
all_spectra[0, :] = fake0
counter = 1
for draw in range(ndraws):
    print("draw #", draw)
    fake1 = np.random.multivariate_normal(np.zeros((N,)), mat)
    for tran in trans:
        ("tran percentage", tran)
        all_spectra[counter, :] = interp_spec(fake0, fake1, tran)
        counter += 1
    fake0 = fake1
        

line, = ax[1].plot(wl23, fake0, "b", drawstyle="steps-mid")

ax[1].set_title("fake residuals")
ax[-1].set_xlabel(r"$\lambda$ [\AA]")
ax[-1].xaxis.set_major_formatter(FSF("%.1f"))
ax[-1].xaxis.set_major_locator(MultipleLocator(0.5))
ax[-1].set_ylim(-0.55, 0.55)
ax[-1].set_xlim(5188.2, 5190)
fig.subplots_adjust(top=0.93, bottom=0.12, right=0.97)

for i,spec in enumerate(all_spectra[1:]):
    line.set_ydata(spec)
    
    fig.canvas.draw_idle()
    fig.savefig("../plots/animation/line/{:0>3d}.png".format(i))









    



[ 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
  0.    0.    0.    0.25  0.5   0.75  1.  ]
shape is (301, 47)
draw # 0
draw # 1
draw # 2
draw # 3
draw # 4
draw # 5
draw # 6
draw # 7
draw # 8
draw # 9
draw # 10
draw # 11
draw # 12
draw # 13
draw # 14



In [53]:

    
#Generating movies using ffmpeg
!ffmpeg -y -r 30 -i ../plots/animation/line/%03d.png -vcodec libx264 -r 30 -pix_fmt yuv420p ../plots/animation/line/line_animation.mp4









    



ffmpeg version 2.2.1 Copyright (c) 2000-2014 the FFmpeg developers
  built on Apr 10 2014 16:26:13 with gcc 4.8.2 (GCC) 20140206 (prerelease)
  configuration: --prefix=/usr --disable-debug --disable-static --enable-avresample --enable-dxva2 --enable-fontconfig --enable-gnutls --enable-gpl --enable-libass --enable-libbluray --enable-libfreetype --enable-libgsm --enable-libmodplug --enable-libmp3lame --enable-libopencore_amrnb --enable-libopencore_amrwb --enable-libopenjpeg --enable-libopus --enable-libpulse --enable-librtmp --enable-libschroedinger --enable-libspeex --enable-libtheora --enable-libv4l2 --enable-libvorbis --enable-libvpx --enable-libx264 --enable-libx265 --enable-libxvid --enable-pic --enable-postproc --enable-runtime-cpudetect --enable-shared --enable-swresample --enable-vdpau --enable-version3 --enable-x11grab
  libavutil      52. 66.100 / 52. 66.100
  libavcodec     55. 52.102 / 55. 52.102
  libavformat    55. 33.100 / 55. 33.100
  libavdevice    55. 10.100 / 55. 10.100
  libavfilter     4.  2.100 /  4.  2.100
  libavresample   1.  2.  0 /  1.  2.  0
  libswscale      2.  5.102 /  2.  5.102
  libswresample   0. 18.100 /  0. 18.100
  libpostproc    52.  3.100 / 52.  3.100
Input #0, image2, from '../plots/animation/line/%03d.png':
  Duration: 00:00:12.00, start: 0.000000, bitrate: N/A
    Stream #0:0: Video: png, rgba, 800x800 [SAR 7874:7874 DAR 1:1], 25 fps, 25 tbr, 25 tbn, 25 tbc
[libx264 @ 0x15abfc0] using SAR=1/1
[libx264 @ 0x15abfc0] using cpu capabilities: MMX2 SSE2Fast SSSE3 SSE4.2
[libx264 @ 0x15abfc0] profile High, level 3.1
[libx264 @ 0x15abfc0] 264 - core 142 r2397 aff928d - H.264/MPEG-4 AVC codec - Copyleft 2003-2014 - http://www.videolan.org/x264.html - options: cabac=1 ref=3 deblock=1:0:0 analyse=0x3:0x113 me=hex subme=7 psy=1 psy_rd=1.00:0.00 mixed_ref=1 me_range=16 chroma_me=1 trellis=1 8x8dct=1 cqm=0 deadzone=21,11 fast_pskip=1 chroma_qp_offset=-2 threads=6 lookahead_threads=1 sliced_threads=0 nr=0 decimate=1 interlaced=0 bluray_compat=0 constrained_intra=0 bframes=3 b_pyramid=2 b_adapt=1 b_bias=0 direct=1 weightb=1 open_gop=0 weightp=2 keyint=250 keyint_min=25 scenecut=40 intra_refresh=0 rc_lookahead=40 rc=crf mbtree=1 crf=23.0 qcomp=0.60 qpmin=0 qpmax=69 qpstep=4 ip_ratio=1.40 aq=1:1.00
Output #0, mp4, to '../plots/animation/line/line_animation.mp4':
  Metadata:
    encoder         : Lavf55.33.100
    Stream #0:0: Video: h264 (libx264) ([33][0][0][0] / 0x0021), yuv420p, 800x800 [SAR 1:1 DAR 1:1], q=-1--1, 15360 tbn, 30 tbc
Stream mapping:
  Stream #0:0 -> #0:0 (png -> libx264)
Press [q] to stop, [?] for help

video:114kB audio:0kB subtitle:0 data:0 global headers:0kB muxing overhead 3.675561%
[libx264 @ 0x15abfc0] frame I:2     Avg QP:17.85  size: 18080
[libx264 @ 0x15abfc0] frame P:96    Avg QP:13.59  size:   617
[libx264 @ 0x15abfc0] frame B:202   Avg QP:14.51  size:   101
[libx264 @ 0x15abfc0] consecutive B-frames: 10.0%  0.7%  0.0% 89.3%
[libx264 @ 0x15abfc0] mb I  I16..4: 36.0% 46.5% 17.6%
[libx264 @ 0x15abfc0] mb P  I16..4:  0.3%  0.3%  0.2%  P16..4:  1.1%  0.4%  0.1%  0.0%  0.0%    skip:97.6%
[libx264 @ 0x15abfc0] mb B  I16..4:  0.0%  0.0%  0.0%  B16..8:  1.0%  0.1%  0.0%  direct: 0.0%  skip:98.8%  L0:52.6% L1:46.7% BI: 0.8%
[libx264 @ 0x15abfc0] 8x8 transform intra:42.5% inter:5.6%
[libx264 @ 0x15abfc0] coded y,uvDC,uvAC intra: 10.0% 14.8% 13.9% inter: 0.1% 0.3% 0.3%
[libx264 @ 0x15abfc0] i16 v,h,dc,p: 64% 34%  2%  0%
[libx264 @ 0x15abfc0] i8 v,h,dc,ddl,ddr,vr,hd,vl,hu: 57%  8% 35%  0%  0%  0%  0%  0%  0%
[libx264 @ 0x15abfc0] i4 v,h,dc,ddl,ddr,vr,hd,vl,hu: 37% 33% 21%  1%  2%  2%  2%  1%  1%
[libx264 @ 0x15abfc0] i8c dc,h,v,p: 83% 11%  6%  0%
[libx264 @ 0x15abfc0] Weighted P-Frames: Y:0.0% UV:0.0%
[libx264 @ 0x15abfc0] ref P L0: 80.1%  1.5% 14.5%  4.0%
[libx264 @ 0x15abfc0] ref B L0: 70.2% 26.7%  3.2%
[libx264 @ 0x15abfc0] ref B L1: 97.8%  2.2%
[libx264 @ 0x15abfc0] kb/s:92.60

M dwarf covariance kernel



In [ ]:

    
def exp_func(x0i, x1i, x0v=None, x1v=None, h=None, amp=None, mu=None, sigma=None):
    x0 = x0v[x0i]
    x1 = x1v[x1i]
    return amp**2/(2 * np.pi * sigma**2) * np.exp(-((x0 - mu)**2 + (x1 - mu)**2)/(2 * sigma**2)) * np.exp(-0.5 * (x0 - x1)**2/h**2)

Matern kernels with a Hanning window function



In [9]:

    
fig = plt.figure()
ax = fig.add_subplot(111)
rs = np.linspace(0,10,num=100)
ax.semilogy(rs, k_3_2(rs, 1), label=r"$k_{3/2}$")
ax.semilogy(rs, hann(rs, 6), label="Hann")
ax.semilogy(rs, k_3_2(rs, 1)* hann(rs, 6), label="comb")
ax.legend()
ax.set_xlim(2, 6)
ax.set_ylim(0, 0.2)
plt.show()

Drawing samples from the covariance matrix



In [5]:

    
len_in = len(xs)
inds = np.arange(len_in)
mat_3_2 = np.fromfunction(k_3_2_hann_func, (len_in,len_in), x0v=xs, x1v=xs, amp=0.8, 
                          l=1.8, dtype=np.int)
y_reals = np.random.multivariate_normal(np.zeros((npoints,)), mat_3_2)
plt.plot(xs, y_reals, "k", lw=0.2)
plt.show()









    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-5-996410b8bdb9> in <module>()
      1 len_in = len(xs)
      2 inds = np.arange(len_in)
----> 3 mat_3_2 = np.fromfunction(k_3_2_hann_func, (len_in,len_in), x0v=xs, x1v=xs, amp=0.8, 
      4                           l=1.8, dtype=np.int)
      5 y_reals = np.random.multivariate_normal(np.zeros((npoints,)), mat_3_2)

NameError: name 'k_3_2_hann_func' is not defined



In [19]:

    
fig,ax = plt.subplots(ncols=2, figsize=(8,4))
ax[0].imshow(mat_3_2, origin="upper", extent=[-20,20,-20,20])
ax[0].set_xlabel(r"$x$")
ax[0].set_title("covariance")
ax[1].plot(xs, y_reals, "k", lw=0.4)
ax[1].set_title("sample")
ax[1].set_xlabel(r"$x$")
ax[1].set_ylabel(r"$y$")
ax[1].set_ylim(-2.5, 2.5)
fig.subplots_adjust(left=0.05, right=0.95, wspace=0.25)
fig.savefig("../plots/matern_matrix.png")
fig.savefig("../plots/matern_matrix.eps")
#plt.show()

Sparse Covariance Matrices



In [6]:

    
#Matern kernel tapered by a Hann window

def sparse_k_3_2(xs, amp, l):
    '''
    Create a sparse matrix using a Matern kernel tapered by a Hann window.
    1. start from the diagonal and keep computing r on the offsets until they are all greater than 6l. 
    2. Then, calculate k_3_2 on the all the diagonals
    3. Initialize a sparse matrix using these diagonals
    '''
    N = len(xs)
    offset = 0
    r0 = 6 * l
    diags = []
    while offset < N:
        #Pairwise calculate rs
        if offset == 0:
            rs = np.zeros_like(xs)
        else:
            rs = np.abs(xs[offset:] - xs[:-offset])
            k = np.empty_like(rs)
        if np.min(rs) >= r0:
            break
        k = (0.5 + 0.5 * np.cos(np.pi * rs/r0)) * amp**2 * (1 + np.sqrt(3) * rs/l) * np.exp(-np.sqrt(3) * rs/l)
        k[rs >= r0] = 0
        diags.append(k)
        offset += 1
        
    #Mirror the diagonals
    front = diags[1:].copy()
    front.reverse()
    diags = front + diags
    offsets = [i for i in range(-offset + 1, offset)]
    return sp.diags(diags, offsets, format="csc")



In [ ]:

    
#To plot residuals from the Matern matrix on top of the data spectrum

model_flux = myOrderModel.get_spectrum()
residuals = fl - model_flux

fig,ax = plt.subplots(nrows=2, figsize=(4,4), sharex=True)
ax[0].plot(wl, fl, label="data")
ax[0].plot(wl, model_flux, "r", label="model")
ax[0].set_ylim(1.e-13, 2.2e-13)
ax[0].legend(loc="lower right")
ax[1].plot(wl, residuals, label="residual")
ax[1].plot(wl_trunc, fl_fake, "k", lw=0.4, label="sample")
ax[1].legend(loc="upper left")
ax[1].set_ylim(-2e-14, 5e-14)
ax[-1].set_xlim(5136, 5139)
ax[-1].set_ylabel(r"$f_\lambda$")
ax[-1].set_xlabel(r"$\lambda$ (\AA)")
fig.subplots_adjust(bottom=0.13, left=0.13, top=0.95, right=0.95)
fig.savefig("../plots/matern_draw.png")
fig.savefig("../plots/matern_draw.svg")
plt.show()



In [144]:

    
ys = np.arange(10)
mysparse = sparse_k_3_2(ys, 1, 1) + sp.eye(len(ys))









    



  (0, 0)	2.0
  (1, 0)	0.450978896606
  (2, 0)	0.104798512644
  (3, 0)	0.0171566215987
  (4, 0)	0.00194193348553
  (5, 0)	0.000112170968055
  (0, 1)	0.450978896606
  (1, 1)	2.0
  (2, 1)	0.450978896606
  (3, 1)	0.104798512644
  (4, 1)	0.0171566215987
  (5, 1)	0.00194193348553
  (6, 1)	0.000112170968055
  (0, 2)	0.104798512644
  (1, 2)	0.450978896606
  (2, 2)	2.0
  (3, 2)	0.450978896606
  (4, 2)	0.104798512644
  (5, 2)	0.0171566215987
  (6, 2)	0.00194193348553
  (7, 2)	0.000112170968055
  (0, 3)	0.0171566215987
  (1, 3)	0.104798512644
  (2, 3)	0.450978896606
  (3, 3)	2.0
  :	:
  (6, 6)	2.0
  (7, 6)	0.450978896606
  (8, 6)	0.104798512644
  (9, 6)	0.0171566215987
  (2, 7)	0.000112170968055
  (3, 7)	0.00194193348553
  (4, 7)	0.0171566215987
  (5, 7)	0.104798512644
  (6, 7)	0.450978896606
  (7, 7)	2.0
  (8, 7)	0.450978896606
  (9, 7)	0.104798512644
  (3, 8)	0.000112170968055
  (4, 8)	0.00194193348553
  (5, 8)	0.0171566215987
  (6, 8)	0.104798512644
  (7, 8)	0.450978896606
  (8, 8)	2.0
  (9, 8)	0.450978896606
  (4, 9)	0.000112170968055
  (5, 9)	0.00194193348553
  (6, 9)	0.0171566215987
  (7, 9)	0.104798512644
  (8, 9)	0.450978896606
  (9, 9)	2.0
  (0, 0)	10.0
  (1, 0)	2.25489448303
  (2, 0)	0.523992563221
  (3, 0)	0.0857831079937
  (4, 0)	0.00970966742763
  (5, 0)	0.000560854840274
  (0, 1)	2.25489448303
  (1, 1)	10.0
  (2, 1)	2.25489448303
  (3, 1)	0.523992563221
  (4, 1)	0.0857831079937
  (5, 1)	0.00970966742763
  (6, 1)	0.000560854840274
  (0, 2)	0.523992563221
  (1, 2)	2.25489448303
  (2, 2)	10.0
  (3, 2)	2.25489448303
  (4, 2)	0.523992563221
  (5, 2)	0.0857831079937
  (6, 2)	0.00970966742763
  (7, 2)	0.000560854840274
  (0, 3)	0.0857831079937
  (1, 3)	0.523992563221
  (2, 3)	2.25489448303
  (3, 3)	10.0
  :	:
  (6, 6)	10.0
  (7, 6)	2.25489448303
  (8, 6)	0.523992563221
  (9, 6)	0.0857831079937
  (2, 7)	0.000560854840274
  (3, 7)	0.00970966742763
  (4, 7)	0.0857831079937
  (5, 7)	0.523992563221
  (6, 7)	2.25489448303
  (7, 7)	10.0
  (8, 7)	2.25489448303
  (9, 7)	0.523992563221
  (3, 8)	0.000560854840274
  (4, 8)	0.00970966742763
  (5, 8)	0.0857831079937
  (6, 8)	0.523992563221
  (7, 8)	2.25489448303
  (8, 8)	10.0
  (9, 8)	2.25489448303
  (4, 9)	0.000560854840274
  (5, 9)	0.00970966742763
  (6, 9)	0.0857831079937
  (7, 9)	0.523992563221
  (8, 9)	2.25489448303
  (9, 9)	10.0



In [6]:

    
def gauss_func(x0i, x1i, x0v=None, x1v=None, amp=None, mu=None, sigma=None):
    x0 = x0v[x0i]
    x1 = x1v[x1i]
    return amp**2/(2 * np.pi * sigma**2) * np.exp(-((x0 - mu)**2 + (x1 - mu)**2)/(2 * sigma**2))



In [9]:

    
def Cfast(xs, amp, mu, sigma, var=1):
    '''Create a sparse covariance matrix using identity and block_diagonal'''
    #In the region of the Gaussian, the matrix will be dense, so just create it as `fromfunction`
    #Above this region, the matrix will be simply Identity.
    
    #The matrix is also symmetric about the diagonal
    
    #Given mu, and the extent of sigma, estimate the data points that are above, in Gaussian, and below
    n_above = np.sum(xs < (mu - 4 * sigma))
    n_below = np.sum(xs > (mu + 4 * sigma))
    ind_in = (xs >= (mu - 4 * sigma)) & (xs <= (mu + 4 * sigma)) #indices to grab the x values 
    len_in = np.sum(ind_in)
    #print(n_above, n_below, len_in)
    #that will be needed to evaluate the Gaussian
    
    if len_in == 0:
        return sp.identity(len(xs), format="csc")
    else:    
        #Create Gaussian, add the sparse diagonal to it
        x_gauss = xs[ind_in]
        gauss_mat = np.fromfunction(gauss_func, (len_in,len_in), x0v=x_gauss, x1v=x_gauss, amp=amp, mu=mu, sigma=sigma, dtype=np.int)
        gauss_mat = gauss_mat + np.identity(len_in)
        
        #plt.imshow(gauss_mat)
        if n_above == 0 and n_below == 0:
            return sp.csc_matrix(gauss_mat)
        elif n_above == 0:
            return sp.block_diag((gauss_mat, sp.identity(n_below)), format="csc")
        elif n_below == 0:
            return sp.block_diag((sp.identity(n_above), gauss_mat), format="csc")
        else:
            return sp.block_diag((sp.identity(n_above), gauss_mat, sp.identity(n_below)), format="csc")

How to initialize a csc_matrix using the list approach



In [8]:

    
data = np.array([3, 3, 3])
ij = np.array([[0, 1, 2],[0, 1 , 2]])
mymat = sp.csc_matrix((data, ij))



In [7]:

    
def k_region(xs, amp, mu, sigma, var=1):
    '''Create a sparse covariance matrix using identity and block_diagonal'''
    #In the region of the Gaussian, the matrix will be dense, so just create it as `fromfunction`
    #and then later turn it into a sparse matrix with size xs x xs
        
    #Given mu, and the extent of sigma, estimate the data points that are above, in Gaussian, and below
    n_above = np.sum(xs < (mu - 4 * sigma))
    n_below = np.sum(xs > (mu + 4 * sigma))
    
    #Create dense matrix and indexes, then convert to lists so that you can pack things in as:
    
    #csc_matrix((data, ij), [shape=(M, N)])
    #where data and ij satisfy the relationship a[ij[0, k], ij[1, k]] = data[k]
    
    len_x = len(xs)
    ind_in = (xs >= (mu - 4 * sigma)) & (xs <= (mu + 4 * sigma)) #indices to grab the x values 
    len_in = np.sum(ind_in)
    #print(n_above, n_below, len_in)
    #that will be needed to evaluate the Gaussian
    
     
    #Create Gaussian matrix fromfunction
    x_gauss = xs[ind_in]
    gauss_mat = np.fromfunction(gauss_func, (len_in,len_in), x0v=x_gauss, x1v=x_gauss, 
                                amp=amp, mu=mu, sigma=sigma, dtype=np.int).flatten()
    
    #Create an index array that matches the Gaussian
    ij = np.indices((len_in, len_in)) + n_above
    ij.shape = (2, -1)
    
    return sp.csc_matrix((gauss_mat, ij), shape=(len_x,len_x))

Create many bad regions



In [13]:

    
S = Cregion(xs, 10, -20, 1) + Cregion(xs, 10, 0, 1) + Cregion(xs, 10, 20, 1) + sp.eye(len(xs))



In [138]:

    
def lnprob_fast(b,m, amp, l):
    A = calcA(b, m)
    S = sparse_k_3_2(xs, amp, l) + sp.eye(len(xs))
    s, logdet = np.linalg.slogdet(S.todense())
    if s <= 0:
        return -np.inf
    lnp =  -0.5 * (A.T.dot(spsolve(S,A)) + logdet)
    return lnp



In [139]:

    
print(lnprob_fast(10, 0.2, 1, 0.1))
print(lnprob_fast(10, 0.2, 1, 1))
print(lnprob_fast(10, 0.2, 1, 10))
print()
print(lnprob_fast(10, 0.2, 10, 0.1))
print(lnprob_fast(10, 0.2, 10, 1))
print(lnprob_fast(10, 0.2, 10, 10))









    



-1883.32528909
-1661.08224326
-1688.58840283

-3224.11425617
-1965.46043112
-1646.16116033

Rather than use a Gaussian for $k(x, x^\prime)$, we can use a more general covariance function, such as the 'squared exponential' but with an added Gaussian taper

$$k(x, x^\prime | h, a, \mu, \sigma) = \exp \left ( \frac{-( x - x^\prime)^2 }{2 h^2} \right ) \frac{a^2}{2 \pi \sigma} \exp \left ( - \frac{[(x - \mu)^2 + (x^\prime - \mu)^2]}{2 \sigma^2}\right )$$

here $h$ is a "bandwidth" that controls the power of the oscillations. If $h$ is small, then there will be high-frequency structure. If $h$ is large, then only low-frequency structure will remain.



In [3]:

    
def exp_func(x0i, x1i, x0v=None, x1v=None, h=None, amp=None, mu=None, sigma=None):
    x0 = x0v[x0i]
    x1 = x1v[x1i]
    return amp**2/(2 * np.pi * sigma**2) * np.exp(-((x0 - mu)**2 + (x1 - mu)**2)/(2 * sigma**2)) * np.exp(-0.5 * (x0 - x1)**2/h**2)