In [1]:
from astropy.io import fits
import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import os
import subprocess
from sklearn.cluster import DBSCAN
from analyzeImage import analyzeImage
import keras
from keras.models import load_model
from __future__ import print_function
%matplotlib inline
%load_ext autoreload
%autoreload 2
In [2]:
gpu_code_path = "../code/gpu/"
real_image_name = "chip_7"
real_image_path = "../../HITS/trimmed_chip_03/Blind15A_02/night_1/" #gpu_code_path+"images/"+real_image_name
psi_image_path = gpu_code_path+"output-images/psi"
phi_image_path = gpu_code_path+"output-images/phi"
In [3]:
def execute(cmd):
popen = subprocess.Popen(cmd, stdout=subprocess.PIPE, universal_newlines=True)
for stdout_line in iter(popen.stdout.readline, ""):
yield stdout_line
popen.stdout.close()
return_code = popen.wait()
if return_code:
raise subprocess.CalledProcessError(return_code, cmd)
In [38]:
paramsFile = open('../code/gpu/debug/parameters.config', 'w')
paramsFile.write(
"""Debug ................ : 1
Image Count .......... : 5
Generate Images ...... : 0
Image Width .......... : 500
Image Height ......... : 500
PSF Sigma ............ : 1.0
Object Brightness .... : 360.0
Object Initial x ..... : 320.0
Object Initial y ..... : 180.0
Velocity x ........... : -1.3
Velocity y ........... : -0.2
Background Level ..... : 1024.0
Background Sigma ..... : 32.0
Mask Threshold ....... : 40.0
Mask Penalty ......... : 15.0
Angles to Search ..... : 80
Velocities to Search . : 20
Minimum Velocity ..... : 300
Maximum Velocity ..... : 450
Write to file ........ : 1
Source Images Path ... : ../../{source}/
Psi Images Path ...... : ../../{psi}/
Phi Images Path....... : ../../{phi}/
""".format( source=real_image_path, psi=psi_image_path, phi=phi_image_path ))
paramsFile.close()
In [5]:
#popen = subprocess.Popen( "./clearImages.sh", stdout=subprocess.PIPE,
# stderr=subprocess.PIPE)
#popen.wait()
#output = popen.stderr.read()
#output += popen.stdout.read()
#print( output)
In [40]:
#for path in execute("./search.sh"):
# print(path, end="")
In [3]:
raw_results = np.genfromtxt('../data/results/HITS1.txt', names=True)
In [4]:
image_mjd = []
for filename in sorted(os.listdir(real_image_path)):
hdulist = fits.open(os.path.join(real_image_path, filename))
image_mjd.append(hdulist[0].header['MJD'])
image_mjd = np.array(image_mjd)
image_times = image_mjd - image_mjd[0]
#image_times*=24.
Load Psi Images
In [5]:
hdulist = fits.open(os.path.join(psi_image_path, os.listdir(psi_image_path)[0]))
num_images = len(os.listdir(psi_image_path))
image_shape = np.shape(hdulist[0].data)
im_psi_array = np.zeros((num_images, image_shape[0], image_shape[1]))
for idx, filename in list(enumerate(sorted(os.listdir(psi_image_path)))):
#print (str('Loaded ' + filename))
image_file = os.path.join(psi_image_path, filename)
hdulist = fits.open(image_file)
im_psi_array[idx] = hdulist[0].data#*mask
Load Phi Images (for potentially making psi/phi stamps)
In [6]:
hdulist = fits.open(os.path.join(phi_image_path, os.listdir(phi_image_path)[0]))
num_images = len(os.listdir(phi_image_path))
image_shape = np.shape(hdulist[0].data)
im_phi_array = np.zeros((num_images, image_shape[0], image_shape[1]))
for idx, filename in list(enumerate(sorted(os.listdir(phi_image_path)))):
# print (str('Loaded ' + filename))
image_file = os.path.join(phi_image_path, filename)
hdulist = fits.open(image_file)
im_phi_array[idx] = hdulist[0].data#*mask
In [7]:
hdulist = fits.open(os.path.join(real_image_path, os.listdir(real_image_path)[0]))
num_images = len(os.listdir(real_image_path))
image_shape = np.shape(hdulist[1].data)
im_array = np.zeros((num_images, image_shape[0], image_shape[1]))
for idx, filename in list(enumerate(sorted(os.listdir(real_image_path)))):
# print( str('Loaded ' + filename))
image_file = os.path.join(real_image_path, filename)
hdulist = fits.open(image_file)
im_array[idx] = hdulist[1].data#*mask
In [8]:
ai = analyzeImage()
In [9]:
model = load_model('../data/kbmod_model.h5')
In [13]:
print(np.shape(raw_results))
results = raw_results[np.where(raw_results['likelihood'] >= 6.0)]
print(np.shape(results))
print(results)
In [14]:
plt.hist(raw_results['likelihood'], bins=np.arange(20))
Out[14]:
In [15]:
filtered_results = ai.filter_results(im_array, results, image_times, model, chunk_size=5000)
In [16]:
print( len(filtered_results) )
filtered_results
Out[16]:
In [17]:
results_to_cluster = filtered_results
arg = dict(eps=0.03, min_samples=1, n_jobs=-1)
clustered_results = ai.clusterResults(results_to_cluster, dbscan_args=arg)#, im_array, image_times)
clustered_results = results_to_cluster[np.array(clustered_results[1], dtype=np.int)]
#best_targets = range(stamp_count)
#best_targets
In [18]:
print( len(clustered_results) )
clustered_results
Out[18]:
In [19]:
#im_psi_array
In [20]:
f_results = clustered_results#filtered_results
rejected = 0
for imNum in range(min(len(f_results),40)):
current = imNum#best_targets[imNum]
fig = plt.figure(figsize=(8,3))
fig.add_subplot(121)
plt.imshow(ai.createPostageStamp(im_array,
list(f_results[['t0_x', 't0_y']][current]),
np.array(list(f_results[['v_x', 'v_y']][current])),
image_times,
[25., 25.])[0],
origin='lower',
#cmap=plt.cm.Greys_r,
interpolation='None')
plt.title(str('#' + str(imNum+1) + ' [x,y] = ' + str(np.round(list(f_results[['t0_x', 't0_y']][current]), 2))
+ ' v = ' + str(np.round(list(f_results[['v_x', 'v_y']][current]),2)) + ' ' +
str(f_results['likelihood'][current])))
fig.add_subplot(122)
ai.plotLightCurves(im_array, f_results[current], image_times)
plt.tight_layout()
traj_coords = ai.calc_traj_coords(f_results[current], image_times)
light_curve = [im_array[x, traj_coords[x,1], traj_coords[x,0]] for x in range(len(im_array))]
if np.max(light_curve) > 50*np.median(light_curve):
rejected+=1
# plt.savefig("stamps/stamp"+str(imNum+1)+".png")
plt.show()
In [62]:
image_times
Out[62]:
In [21]:
rejected
Out[21]:
In [ ]: