In [1]:
%matplotlib inline
%config InlineBackend.figure_format = "retina"
from matplotlib import rcParams
rcParams["savefig.dpi"] = 100
rcParams["font.size"] = 20
In [171]:
import kplr
import numpy as np
import matplotlib.pyplot as pl
from simplexy import simplexy
In [113]:
client = kplr.API()
tpf = client.k2_star(202137899).get_target_pixel_files()[0]
In [277]:
data = tpf.read()
time = data["TIME"]
flux = data["FLUX"]
ferr = data["FLUX_ERR"]
quality = data["QUALITY"]
m = np.all((flux == 0.0) + np.isnan(flux), axis=(1, 2)) + np.isnan(time) + (quality > 0)
t = time[~m]
i = np.argmax(np.diff(t))
i = np.arange(len(time))[~m][i+1]
i = -1000
time = time[i:]
flux = flux[i:]
ferr = ferr[i:]
quality = quality[i:]
print(i, time[0])
In [278]:
def run_one(img):
m = np.isfinite(img)
m[m] *= (img[m] > 0.0)
if not np.any(m):
return None
img[~m] = np.median(img[m])
return np.array([list(row) for row in simplexy(img, )])
Generate a graph of coordinate locations.
In [279]:
graph = []
for i in range(len(flux)):
if quality[i] != 0:
graph.append([])
continue
img = flux[i]
r = run_one(img)
graph.append(r if r is not None else [])
Traverse the graph greedily.
In [289]:
nt = len(graph)
ns = min(map(len, filter(len, graph)))
coords = np.nan + np.zeros((nt, ns, 4))
metric = np.array([1.0, 1.0, 1e-8])
current = None
for t, node in enumerate(graph):
if not len(node):
continue
if current is None:
current = node[:ns, :]
coords[t] = current
continue
# Compute the set of distances between this node and the current position.
r = np.sum((node[:, None, :3] - current[None, :, :3])**2 * metric[None, None, :], axis=-1)
# Loop over the permutations and greedily choose the best update.
rows, cols = np.arange(r.shape[0]), np.arange(r.shape[1])
update = np.nan + np.zeros(ns)
while any(np.isnan(update)):
row, col = np.unravel_index(np.argmin(r), r.shape)
update[cols[col]] = rows[row]
r = np.delete(r, row, axis=0)
r = np.delete(r, col, axis=1)
rows = np.delete(rows, row)
cols = np.delete(cols, col)
update = np.array(update, dtype=int)
# Compute the total cost.
cost = np.sum((current[:, :3] - node[update, :3])**2 * metric[None, :])
if cost > 10.0:
continue
# Update the current locations.
current = node[update]
coords[t] = current
Generate a mask of times that are OK.
In [294]:
good_times = np.all(np.isfinite(coords), axis=(1, 2))
Find the mean motion of the frame.
In [295]:
center = np.array(coords[:, 0, :2])
center -= np.mean(center[good_times], axis=0)
# center = np.sum(coords[:, :, :2] * coords[:, :, 2][:, :, None], axis=1) / np.sum(coords[:, :, 2], axis=1)[:, None]
dx = coords[:, :, :2] - center[:, None, :]
offsets = np.median(dx[np.all(np.isfinite(dx), axis=(1, 2))], axis=0)
In [296]:
pl.plot(time, coords[:, 1, 2], ".k")
Out[296]:
In [299]:
pl.scatter(center[:, 0], center[:, 1], c=time)
Out[299]:
Make a movie!
In [301]:
i = 0
for t, img in enumerate(flux):
m = np.isfinite(img)
m[m] *= (img[m] > 0.0)
if not np.any(m):
continue
img[~m] = np.median(img[m])
pl.clf()
pl.imshow(np.log(img.T), cmap="gray_r", interpolation="nearest")
pl.plot(coords[t, :, 0], coords[t, :, 1], "+b")
pl.plot(center[t, 0] + offsets[:, 0], center[t, 1] + offsets[:, 1], "+r")
pl.title("EPIC 202137899")
pl.savefig("frames/{0:05d}.png".format(i))
i += 1
In [302]:
offsets
Out[302]:
In [ ]: