In [1]:
%matplotlib inline
from matplotlib import rcParams
rcParams["savefig.dpi"] = 150
In [191]:
import h5py
import json
import fitsio
import numpy as np
import pandas as pd
from astropy.wcs import WCS
from itertools import product
import matplotlib.pyplot as pl
from scipy.linalg import cho_factor, cho_solve
from matplotlib.colors import LinearSegmentedColormap
from sklearn.decomposition import FastICA, RandomizedPCA
In [198]:
with open("k2.wcs.json", "r") as f:
wcs = WCS(json.load(f))
In [199]:
df = pd.read_csv("K2_var.dat", delim_whitespace=True, names=[
"ID", "RA", "DEC", "X", "Y", "Mag_Sch", "Mag_K2", "Phot_K2", "Clas_K2", "Period",
"ID_Nar", "ID_Hu", "ID_Hu", "ID_Kim", "ID_Mei", "ID_Moc", "B", "V", "R",
"J_2MASS", "H_2MASS", "K_2MASS",
])
df = df[df.Clas_K2 == 1]
In [228]:
ra = np.array(df.RA)
dec = np.array(df.DEC)
x, y = wcs.all_world2pix(ra, dec, 1)
df["k2_cx"] = x
df["k2_cy"] = y
print(x.min(), x.max(), y.min(), y.max())
In [469]:
df[df.J_2MASS > 14].sort("J_2MASS")
Out[469]:
In [195]:
with h5py.File("data/k2.h5", "r") as f:
time = f["time"][:]
i = np.arange(len(time))[np.isnan(time)][-1] + 100
time = time[i+1:]
data = f["frames"][i+1:, :, :]
data[data == 0] = np.nan
In [196]:
m = np.isfinite(time) & np.any(np.isfinite(data), axis=(1, 2))
time = time[m]
data = data[m]
In [206]:
pl.imshow(np.log(data[-1]), cmap="gray", interpolation="nearest")
pl.plot(x, y, ".r", ms=2);
In [67]:
scaled = data # / np.median(data, axis=0)[None, :, :]
In [68]:
data.shape
Out[68]:
In [538]:
X, Y = np.meshgrid(range(scaled.shape[1]), range(scaled.shape[2]), indexing="ij")
m = np.any(np.isfinite(scaled), axis=0)
pred = np.zeros_like(scaled)
count = np.zeros_like(scaled[0], dtype=int)
hw = 5
ex = hw + 10
cen = [258.794798, 272.821687]
xmn, xmx = int(cen[1] - 5), int(cen[1] + 5)
ymn, ymx = int(cen[0] - 4), int(cen[0] + 4)
# xmn, xmx = 575, 585
# ymn, ymx = 515, 525
# for i, j in product(range(0, scaled.shape[1], hw+1), range(0, scaled.shape[1], hw+1)):
for i, j in product(range(xmn, xmx, hw+1), range(ymn, ymx, hw+1)):
# Select target and predictor pixels.
window = m & (i - hw <= X) & (X <= i + hw) & (j - hw <= Y) & (Y <= j + hw)
r2 = (X-i)**2 + (Y-j)**2
others = m & ((i - ex > X) | (X > i + ex) | (j - ex > Y) | (Y > j + ex)) & (r2 < 25**2)
print(i, j, window.sum(), others.sum())
if not (np.any(window) and np.any(others)):
continue
A = np.concatenate((scaled[:, others], np.ones((len(scaled), 1))), axis=1)
ATA = np.dot(A.T, A)
ATA[np.diag_indices_from(ATA)] += 1e-2
factor = cho_factor(ATA, overwrite_a=True)
# y = np.concatenate((scaled[:, window], np.zeros((len(emp), window.sum()))), axis=0)
y = scaled[:, window]
w = cho_solve(factor, np.dot(A.T, y))
pred[:, window] += np.dot(A[:len(scaled)], w)
count[window] += 1
In [539]:
pred[:, count > 0] /= count[count > 0]
In [540]:
d = (scaled[:, xmn:xmx, ymn:ymx] - pred[:, xmn:xmx, ymn:ymx])
np.unravel_index(np.argmax(np.abs(d)), d.shape)
Out[540]:
In [541]:
n = 100
pl.imshow(data[n, xmn:xmx, ymn:ymx] - pred[n, xmn:xmx, ymn:ymx], cmap="gray", interpolation="nearest")
pl.colorbar()
Out[541]:
In [542]:
d = data[:, xmn:xmx, ymn:ymx] - pred[:, xmn:xmx, ymn:ymx]
block = d.reshape((len(d), -1))
# block = block / np.median(block, axis=0) - 1
In [543]:
pl.plot(time, np.sum(block, axis=1), ".k")
# pl.plot(time, np.sum(block, axis=1), ".k")
Out[543]:
In [550]:
model = FastICA(n_components=2)
model.fit(block.T)
Out[550]:
In [551]:
n = len(model.components_)
print(n)
In [552]:
pl.figure(figsize=(5, 10))
pl.plot(model.components_.T + 1e-4*np.arange(n));
In [553]:
fig, axes = pl.subplots(1, 2, figsize=(8, 4))
v = model.components_[1]
comp = np.abs((np.dot(v, block)).reshape(d[0].shape))
axes[0].imshow(np.log(data[0, xmn:xmx, ymn:ymx]), cmap="gray_r", interpolation="nearest")
axes[1].imshow(comp, cmap="gray_r", interpolation="nearest")
Out[553]:
In [557]:
c = np.dot(v, block)
c /= np.sum(c)
t = time % 5.862149
# y = np.sum(block, axis=1)
# pl.plot(t, y / np.median(np.abs(np.diff(y))), ".g")
y = np.sum(c * block, axis=1)
pl.plot(t, y / np.median(np.abs(np.diff(y))), ".k")
Out[557]:
In [556]:
pl.figure()
pl.imshow(d[0], cmap="gray_r", interpolation="nearest")
pl.figure()
pl.imshow(d[0], cmap="gray_r", interpolation="nearest")
pl.imshow(comp, cmap="Blues", interpolation="nearest", alpha=0.5);
In [ ]:
In [ ]: