In [1]:
import pandas as pd
import numpy as np
import os,sys
import random
import re
import seaborn as sns
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.colors import BoundaryNorm
import collections
from keras import models
from keras import layers
from keras.optimizers import SGD, Adadelta
from keras.callbacks import EarlyStopping, ModelCheckpoint
import datetime
import h5py
from sklearn.metrics import roc_curve, auc
from mpl_toolkits.axes_grid.inset_locator import inset_axes, InsetPosition
%matplotlib inline
In [2]:
def show_images(indice, ref_images, new_images, sub_images, title=None, Npanel=None, labels=None, classes=None, labelAsValue=True, labelAsFloat=False, save_name=None):
if Npanel == None:
Npanel = (len(indice)-1)/10 + 1
for ii in range(Npanel):
fig = plt.figure(figsize=(14, 14))
if title != None:
fig.suptitle(title+' (%d)' % (ii+1), fontsize=16)
n = 0
for i in indice[ii*10:]:
ref_ax = fig.add_subplot(5, 6, 3*n+0+1)
ref_ax.set_xticks([])
ref_ax.set_yticks([])
ref_ax.imshow(ref_images[i].reshape(29,29),
cmap=plt.cm.gray,
#norm=BoundaryNorm(np.linspace(ref_images[i].min(), ref_images[i].min()+0.2, 256), ncolors=256),
interpolation='nearest')
new_ax = fig.add_subplot(5, 6, 3*n+1+1)
new_ax.set_xticks([])
new_ax.set_yticks([])
new_ax.imshow(new_images[i].reshape(29,29),
cmap=plt.cm.gray,
#norm=BoundaryNorm(np.linspace(ref_images[i].min(), ref_images[i].mean()+500, 256), ncolors=256),
interpolation='nearest')
sub_ax = fig.add_subplot(5, 6, 3*n+2+1)
sub_ax.set_xticks([])
sub_ax.set_yticks([])
sub_ax.imshow(sub_images[i].reshape(29,29),
cmap=plt.cm.gray,
#norm=BoundaryNorm(np.linspace(ref_images[i].min(), ref_images[i].min()+2.0, 256), ncolors=256),
interpolation='nearest')
if labels is not None and classes is not None:
ref_ax.set_title('ref:{:1.0f}/{:.2f} ({})'.format(labels[i], classes[i], i))
new_ax.set_title('new:{:1.0f}/{:.2f} ({})'.format(labels[i], classes[i], i))
sub_ax.set_title('sub:{:1.0f}/{:.2f} ({})'.format(labels[i], classes[i], i))
#ax.set_title('%d' % (i))
elif labels is not None:
if labelAsValue:
if labelAsFloat:
ref_ax.set_title('ref:%f' % (labels[i]))
new_ax.set_title('new:%f' % (labels[i]))
sub_ax.set_title('sub:%f' % (labels[i]))
else:
ref_ax.set_title('ref:{}'.format(labels[i]))
new_ax.set_title('new:{}'.format(labels[i]))
sub_ax.set_title('sub:{}'.format(labels[i]))
else:
ref_ax.set_title('ref:%d' % (i))
new_ax.set_title('new:%d' % (i))
sub_ax.set_title('sub:%d' % (i))
else:
ref_ax.set_title('ref:%d' % (i))
new_ax.set_title('new:%d' % (i))
sub_ax.set_title('sub:%d' % (i))
n += 1
if n == 10:
break
if save_name is not None: plt.savefig(save_name+'_{}.png'.format(ii))
plt.show()
plt.close()
In [3]:
def norm_image(ref_image_nd, new_image_nd, sub_image_nd):
ref_image_nd = ref_image_nd.reshape(ref_image_nd.shape[0],ref_image_nd.shape[1])
new_image_nd = new_image_nd.reshape(new_image_nd.shape[0],new_image_nd.shape[1])
sub_image_nd = sub_image_nd.reshape(sub_image_nd.shape[0],sub_image_nd.shape[1])
ref_image_nd = (ref_image_nd-ref_image_nd.min())/(ref_image_nd.max()-ref_image_nd.min())
new_image_nd = (new_image_nd-new_image_nd.min())/(new_image_nd.max()-new_image_nd.min())
sub_image_nd = (sub_image_nd-sub_image_nd.min())/(sub_image_nd.max()-sub_image_nd.min())
return [ref_image_nd, new_image_nd, sub_image_nd]
Use Model ver.5 (New)¶
In [4]:
root_directory = '/home/tomoesn/script/CNN/cnn_model_v5'
model_dist = {}
det_q_list = [1,2,3,4]
det_num_list = [11,12,13,14,15,16,
21,22,23,24,25,26,
31,32,33,34,35,
41,42,43,44]
load_det_num = [q*100+num for q in det_q_list for num in det_num_list]
for det_num in load_det_num:
print 'Load {} model.'.format(det_num)
model_file_path = os.path.join(root_directory, str(det_num), 'model{}.hdf5'.format(det_num))
model_dist[det_num] = models.load_model(model_file_path)
In [5]:
csv_path = '/home/hamasaki/hamasaki-2/TestTransient/release/data/real_transient_2019-09-01_2019-11-01.csv'
npy_path = '/home/hamasaki/hamasaki-2/TestTransient/release/data/real_transient_2019-09-01_2019-11-01.npy'
real_params = pd.read_csv(csv_path)
real_images = np.load(npy_path)
print real_params.shape
print real_images.shape
In [6]:
pd.set_option('display.max_columns',80)
In [7]:
real_params[:6]
Out[7]:
In [8]:
index = [0,1,2,3,4,5]
images_ch_list = np.split(real_images, 3, axis=3)
print index
show_images(index,*images_ch_list, title='TNS Transient')
In [9]:
index = np.array(real_params[real_params['ps_match']==False].index)
images_ch_list = np.split(real_images, 3, axis=3)
print index.shape
print index
show_images(index[10:20],*images_ch_list, title='TNS Transient without star')
Let's look at REAL (TNS) Transient.¶
In [10]:
real_norm_images = np.array([np.stack(norm_image(*np.split(images, 3, axis=2)), axis=-1) for images in real_images])
In [11]:
real_prob = []
real_class = []
for det_num, images in zip(real_params['tom_det_id'].values.astype('int64'), real_norm_images):
images = images.reshape(1, images.shape[0], images.shape[1], images.shape[2])
predict_prob = model_dist[det_num].predict_proba(images, verbose=False)
predict_class = model_dist[det_num].predict_classes(images, verbose=False)
real_prob.append(predict_prob[0,0])
real_class.append(predict_class[0,0])
real_prob = np.array(real_prob)
real_class = np.array(real_class)
In [12]:
real_params['ana_cnn_cand'] = real_class
real_params['ana_cnn_prob'] = real_prob
In [13]:
real_params[:5]
Out[13]:
In [14]:
detected_real_params = real_params[(real_params['ana_cnn_cand']==1)&(real_params['ps_match']==False)&
(real_params['sub_value_zero']==False)&(real_params['sub_mask']==False)&
(real_params['sub_edge']==False)&(real_params['sub_saturation']==False)]
In [15]:
print real_params[real_params['ps_match']==False].shape
print detected_real_params.shape
In [16]:
not_detected_real_params = real_params[(real_params['ana_cnn_cand']==0)&(real_params['ps_match']==False)&
(real_params['sub_value_zero']==False)&(real_params['sub_mask']==False)&
(real_params['sub_edge']==False)&(real_params['sub_saturation']==False)]
In [17]:
index = np.array(detected_real_params.index)
images_ch_list = np.split(real_images, 3, axis=3)
print index
show_images(index[:10],*images_ch_list, title='Detected')
In [18]:
index = np.array(not_detected_real_params.index)
images_ch_list = np.split(real_images, 3, axis=3)
print index
show_images(index,*images_ch_list, title='Not Detected')
In [19]:
real_params[real_params['ps_type']=='star']
Out[19]:
In [20]:
real_params[real_params['ps_id']==0]
Out[20]:
In [21]:
print real_params[real_params['tns_name'].str.contains('SN')].shape
print real_params.shape
In [22]:
star_params = real_params[(real_params['ps_type']=='star')&(real_params['ps_id']!=0)]
ext_params = real_params[(real_params['ps_type']=='extended')&(real_params['ps_id']!=0)]
In [23]:
print star_params.shape
print ext_params.shape
In [24]:
fig = plt.figure(figsize=(15,7))
fig.suptitle('PS_dist vs PS_rmag 2019/09/01-2019/11/01',fontsize=18)
ax = fig.add_subplot(121)
ax.plot(star_params['ps_dist'].values*3600, star_params['ps_rmag'].values,
'.', label='TNS (AT)', color='tab:blue')
ax.plot(star_params[star_params['tns_name'].str.contains('SN')]['ps_dist'].values*3600,
star_params[star_params['tns_name'].str.contains('SN')]['ps_rmag'].values,
'x', label='TNS (SN)', color='tab:red')
ax.vlines([3],20,16,'gray',linestyles='dashed',label='threshold')
ax.vlines([10],10,16,'gray',linestyles='dashed')
ax.hlines([16],3,10,'gray',linestyles='dashed')
ax.text(1,11,'PSstar==1', size=20, color='green')
ax.set_xlim(0, 60)
ax.set_ylim(20, 10)
#ax.invert_yaxis()
ax.set_xlabel('distance (arcsec)',fontsize=15)
ax.set_ylabel('r_mag',fontsize=15)
ax.tick_params(labelsize=12)
plt.legend(fontsize=12)
ax = fig.add_subplot(122)
ax.plot(ext_params['ps_dist'].values*3600, ext_params['ps_rmag'].values,
'.', label='TNS (AT)', color='tab:blue')
ax.plot(ext_params[ext_params['tns_name'].str.contains('SN')]['ps_dist'].values*3600,
ext_params[ext_params['tns_name'].str.contains('SN')]['ps_rmag'].values,
'x', label='TNS (SN)', color='tab:red')
ax.vlines([3],20,16,'gray',linestyles='dashed',label='threshold')
ax.vlines([10],10,16,'gray',linestyles='dashed')
ax.hlines([16],3,10,'gray',linestyles='dashed')
ax.set_xlim(0, 60)
ax.set_ylim(20, 10)
#ax.invert_yaxis()
ax.set_xlabel('distance (arcsec)',fontsize=15)
ax.set_ylabel('r_mag',fontsize=15)
ax.tick_params(labelsize=12)
ax.text(1,11,'EXT==1', size=20, color='green')
plt.legend(fontsize=12)
plt.savefig('ps_dist_vs_ps_rmag.png',dpi=150)
plt.show()
plt.close()
In [25]:
index = np.array(real_params[real_params['ps_id']==0].index)
print index
show_images(index,*np.split(real_images, 3, axis=3), title='PS Over 60 arcsec')
In [26]:
index = star_params[(star_params['ps_dist']<(3.0/3600))&(star_params['tns_name'].str.contains('AT'))].index.values
param=star_params[(star_params['ps_dist']<(3.0/3600))&(star_params['tns_name'].str.contains('AT'))]
print index
show_images(index[:10],*np.split(real_images, 3, axis=3), labelAsFloat=True,
labels=param['ana_cnn_prob'],
#save_name='psstar_and_dist_less_than_3arcsec',
title='AT\nPSstar==1 & Dist<3arcsec')
In [27]:
print param.shape
print param['ana_cnn_cand'].sum()
In [28]:
index = ext_params[(ext_params['ps_dist']<(3.0/3600))&(ext_params['tns_name'].str.contains('AT'))].index.values
param=ext_params[(ext_params['ps_dist']<(3.0/3600))&(ext_params['tns_name'].str.contains('AT'))]
print index
show_images(index[:10],*np.split(real_images, 3, axis=3), labelAsFloat=True,
labels=param['ana_cnn_prob'],
#save_name='psext_and_dist_less_than_3arcsec',
title='AT\nEXT==1 & Dist<3arcsec')
Let's look at BOGUS Transient.¶
In [29]:
csv_path = '/home/hamasaki/hamasaki-2/TestTransient/release/data/bogus_transient_2020-02-01_2020-03-01.csv'
npy_path = '/home/hamasaki/hamasaki-2/TestTransient/release/data/bogus_transient_2020-02-01_2020-03-01.npy'
bogus_params = pd.read_csv(csv_path)
bogus_images = np.load(npy_path)
print bogus_params.shape
print bogus_images.shape
In [30]:
import psycopg2
In [31]:
def inquire_tns(connection, bogus_params, temp_table_name='bogus_transient_tns_match'):
ra = bogus_params['sub_ra'].values
dec = bogus_params['sub_dec'].values
mjd = bogus_params['tom_mjd'].values
ra_dec_mjd = np.stack([ra,dec,mjd], axis=-1)
ra_dec_mjd = tuple(map(tuple,ra_dec_mjd))
ra_dec_mjd = ','.join(map(str, ra_dec_mjd))
with connection.cursor() as cur:
query = 'CREATE TEMPORARY TABLE {} (ra DOUBLE PRECISION, dec DOUBLE PRECISION, mjd DOUBLE PRECISION)'.format(temp_table_name)
cur.execute(query)
query = 'INSERT INTO {} (ra,dec,mjd) VALUES {}'.format(temp_table_name, ra_dec_mjd)
cur.execute(query)
connection.commit()
with connection.cursor('match_tns_table') as cur: # ---!!!CAUTION!!!--- Heavy Processing
cur.itersize = 10000
# Search tns
min_date = '2019-11-01'
max_date = '2020-02-01'
#query = 'SELECT ss.* FROM {} AS t LEFT JOIN LATERAL (SELECT s."tnsId", s."ra", s."dec", s."discMag", s."name", s."discDate", q3c_dist(t.ra, t.dec, s.ra, s.dec) FROM tns AS s WHERE s."discDate" > \'{}\' AND s."discDate" < \'{}\' AND q3c_join(t.ra, t.dec, s.ra, s.dec, 1./60) ORDER BY q3c_dist(t.ra, t.dec, s.ra, s.dec) ASC LIMIT 1 ) AS ss ON true'.format(temp_table_name, min_date, max_date)
query = 'SELECT ss.* FROM {} AS t LEFT JOIN LATERAL (SELECT s."tnsId", s."ra", s."dec", s."discMag", s."name", s."discDate", q3c_dist(t.ra, t.dec, s.ra, s.dec) FROM tns AS s WHERE q3c_join(t.ra, t.dec, s.ra, s.dec, 60.0/3600.0) ORDER BY q3c_dist(t.ra, t.dec, s.ra, s.dec) ASC LIMIT 1 ) AS ss ON true'.format(temp_table_name)
cur.execute(query)
responses = cur.fetchall()
data = [[0 if num is None else num for num in ps_param] for ps_param in responses]
tns_df = pd.DataFrame(data, columns=['tns_id','tns_ra','tns_dec','tns_disc_mag','tns_name','tns_disc_date','tns_dist'])
return tns_df
In [32]:
psql_login = 'host=gwdata dbname=supernova user=tomoesn password='
with psycopg2.connect(psql_login) as connection:
tns_df = inquire_tns(connection, bogus_params)
In [33]:
tns_df[tns_df['tns_id']!=0][:5]
Out[33]:
In [34]:
bogus_params[:3]
Out[34]:
In [35]:
index = [0,1,2,3,4,5]
print index
show_images(index,*np.split(bogus_images, 3, axis=3), title='BOGUS Transient')
In [36]:
fig = plt.figure(figsize=(7,7))
fig.suptitle('BOGUS Coordinate Distribution',fontsize=18)
ax = fig.add_subplot(111)
ax.plot(bogus_params['sub_ra'].values, bogus_params['sub_dec'].values,
'.', label='BOGUS Transient', color='tab:blue')
ax.set_xlim(0, 360)
ax.set_ylim(-20, 90)
#ax.invert_yaxis()
ax.set_xlabel('ra (deg)',fontsize=15)
ax.set_ylabel('dec (deg)',fontsize=15)
ax.tick_params(labelsize=12)
plt.legend(fontsize=12)
#plt.savefig('ps_dist_vs_ps_rmag.png',dpi=150)
plt.show()
plt.close()
In [37]:
bogus_norm_images = np.array([np.stack(norm_image(*np.split(images, 3, axis=2)), axis=-1) for images in bogus_images])
In [38]:
bogus_prob = []
bogus_class = []
for det_num, images in zip(bogus_params['tom_det_id'].values.astype('int64'), bogus_norm_images):
images = images.reshape(1, images.shape[0], images.shape[1], images.shape[2])
predict_prob = model_dist[det_num].predict_proba(images, verbose=False)
predict_class = model_dist[det_num].predict_classes(images, verbose=False)
bogus_prob.append(predict_prob[0,0])
bogus_class.append(predict_class[0,0])
bogus_prob = np.array(bogus_prob)
bogus_class = np.array(bogus_class)
bogus_params['ana_cnn_cand'] = bogus_class
bogus_params['ana_cnn_prob'] = bogus_prob
Let's look at REAL&BOGUS Transient.¶
In [39]:
def check_param_cand(param):
param_cand_flag = True
if param['sub_value_zero'] == True : param_cand_flag = False
if param['sub_saturation'] == True : param_cand_flag = False
if param['sub_mask'] == True : param_cand_flag = False
if param['sub_edge'] == True : param_cand_flag = False
return param_cand_flag
In [40]:
real_params['ana_param_cand'] = real_params.apply(check_param_cand, axis=1)
bogus_params['ana_param_cand'] = bogus_params.apply(check_param_cand, axis=1)
In [41]:
def check_ps_type(param):
if param['ps_id'] == 0: return 'none'
return param['ps_type']
In [42]:
real_params['ps_type'] = real_params.apply(check_ps_type, axis=1)
bogus_params['ps_type'] = bogus_params.apply(check_ps_type, axis=1)
In [43]:
bogus_params[bogus_params['ana_param_cand']==True][:5]
Out[43]:
In [44]:
real_params[:5]
Out[44]:
In [45]:
print bogus_params[bogus_params['ana_param_cand']==True].shape
print real_params[real_params['ana_param_cand']==True].shape
In [46]:
label = np.concatenate([np.ones(real_params[real_params['ana_param_cand']==True]['ana_cnn_prob'].shape[0]),
np.zeros(bogus_params[bogus_params['ana_param_cand']==True]['ana_cnn_prob'].shape[0])],
axis=0)
prob = np.concatenate([real_params[real_params['ana_param_cand']==True]['ana_cnn_prob'].values,
bogus_params[bogus_params['ana_param_cand']==True]['ana_cnn_prob'].values],
axis=0)
fpr, tpr, thresholds = roc_curve(label, prob)
val_auc = auc(fpr, tpr)
fig = plt.figure(figsize=(8,6))
fig.suptitle("ROC curve for Test Set (AUC=%6.4f)" % (val_auc), fontsize=16)
ax = fig.add_subplot(111)
ax.plot(fpr, tpr, color='tab:red', label='AUC={:>6.3f}'.format(val_auc))
ax.set_xlabel("False Positive Rate",fontsize=15)
ax.set_ylabel("True Positive Rate",fontsize=15)
ax.set_ylim(-0.05,1.05)
ax.set_xlim(-0.05,1.05)
ax.tick_params(labelsize=12)
#ax.set_xscale('log')
ax.grid()
ax.set_aspect('equal')
iaxes = inset_axes(ax, width="50%", height="50%", loc=4)
ip = InsetPosition(ax, [0.3,0.1,0.6,0.6])
iaxes.set_axes_locator(ip)
iaxes.plot(fpr, tpr, color='tab:red')
iaxes.set_xticks(np.arange(0.0, 0.2, 0.050))
iaxes.set_xticks(np.arange(0.0, 0.2, 0.010), minor=True)
iaxes.set_yticks(np.arange(0.70, 1.0, 0.020))
iaxes.set_yticks(np.arange(0.70, 1.0, 0.010), minor=True)
iaxes.set_xlim(0.00, 0.15)
iaxes.set_ylim(0.80, 0.98)
iaxes.tick_params(labelsize=12)
iaxes.grid()
iaxes.set_aspect('equal')
plt.savefig('new_model_roc_curve.png')
plt.show()
plt.close()
Exclude 'variable star'¶
In [47]:
real_params[~((real_params['ps_type']=='star')&(real_params['ps_match']==True))&(real_params['ana_param_cand']==True)][:5]
Out[47]:
In [48]:
bogus_params[~((bogus_params['ps_type']=='star')&(bogus_params['ps_match']==True))&(bogus_params['ana_param_cand']==True)][:5]
Out[48]:
In [49]:
label = np.concatenate([np.ones(real_params[~((real_params['ps_type']=='star')&(real_params['ps_match']==True))&(real_params['ana_param_cand']==True)]['ana_cnn_prob'].shape[0]),
np.zeros(bogus_params[~((bogus_params['ps_type']=='star')&(bogus_params['ps_match']==True))&(bogus_params['ana_param_cand']==True)]['ana_cnn_prob'].shape[0])],
axis=0)
prob = np.concatenate([real_params[~((real_params['ps_type']=='star')&(real_params['ps_match']==True))&(real_params['ana_param_cand']==True)]['ana_cnn_prob'].values,
bogus_params[~((bogus_params['ps_type']=='star')&(bogus_params['ps_match']==True))&(bogus_params['ana_param_cand']==True)]['ana_cnn_prob'].values],
axis=0)
fpr, tpr, thresholds = roc_curve(label, prob)
val_auc = auc(fpr, tpr)
fig = plt.figure(figsize=(8,6))
fig.suptitle("ROC curve for Test Set (AUC=%6.4f)" % (val_auc), fontsize=16)
ax = fig.add_subplot(111)
ax.plot(fpr, tpr, color='tab:red', label='AUC={:>6.3f}'.format(val_auc))
ax.set_xlabel("False Positive Rate",fontsize=15)
ax.set_ylabel("True Positive Rate",fontsize=15)
ax.set_ylim(-0.05,1.05)
ax.set_xlim(-0.05,1.05)
ax.tick_params(labelsize=12)
#ax.set_xscale('log')
ax.grid()
ax.set_aspect('equal')
iaxes = inset_axes(ax, width="50%", height="50%", loc=4)
ip = InsetPosition(ax, [0.3,0.1,0.6,0.6])
iaxes.set_axes_locator(ip)
iaxes.plot(fpr, tpr, color='tab:red')
iaxes.set_xticks(np.arange(0.0, 0.2, 0.050))
iaxes.set_xticks(np.arange(0.0, 0.2, 0.010), minor=True)
iaxes.set_yticks(np.arange(0.70, 1.0, 0.020))
iaxes.set_yticks(np.arange(0.70, 1.0, 0.010), minor=True)
iaxes.set_xlim(0.00, 0.15)
iaxes.set_ylim(0.80, 0.98)
iaxes.tick_params(labelsize=12)
iaxes.grid()
iaxes.set_aspect('equal')
plt.savefig('new_model_roc_curve.png')
plt.show()
plt.close()
In [50]:
real = real_params[~((real_params['ps_type']=='star')&(real_params['ps_match']==True))&
(real_params['ana_param_cand']==True)]
bogus = bogus_params[~((bogus_params['ps_type']=='star')&(bogus_params['ps_match']==True))&
(bogus_params['ana_param_cand']==True)]
tp = real_params[~((real_params['ps_type']=='star')&(real_params['ps_match']==True))&
(real_params['ana_param_cand']==True)&
(real_params['ana_cnn_prob']>=0.5)]
tn = bogus_params[~((bogus_params['ps_type']=='star')&(bogus_params['ps_match']==True))&
(bogus_params['ana_param_cand']==True)&
(bogus_params['ana_cnn_prob']<0.5)]
fp = bogus_params[~((bogus_params['ps_type']=='star')&(bogus_params['ps_match']==True))&
(bogus_params['ana_param_cand']==True)&
(bogus_params['ana_cnn_prob']>=0.5)]
print 'True Positive : {:>5}/{:>5} = {:>5.3f}'.format(tp.shape[0], real.shape[0], float(tp.shape[0])/real.shape[0])
print 'False Positive : {:>5}/{:>5} = {:>5.3f}'.format(fp.shape[0], bogus.shape[0], float(fp.shape[0])/bogus.shape[0])
In [53]:
fig = plt.figure(figsize=(12,6))
fig.suptitle("TPR, FPR vs. Threshold", fontsize=18)
ax = fig.add_subplot(121)
ax.plot(thresholds, fpr, color='tab:red')
ax.set_xlabel("Threshold", fontsize=15)
ax.set_ylabel("False Positive Rate", fontsize=15)
ax.set_ylim(-0.05,1.05)
ax.set_xlim(-0.05,1.05)
ax.tick_params(labelsize=12)
ax.set_aspect('equal')
iaxes = inset_axes(ax, width="50%", height="50%", loc=4)
ip = InsetPosition(ax, [0.20,0.65,0.7,0.3])
iaxes.set_axes_locator(ip)
iaxes.tick_params(labelsize=12)
iaxes.plot(thresholds, fpr, color='tab:red')
iaxes.set_xticks(np.arange(0.0, 1.1, 0.100))
iaxes.set_xticks(np.arange(0.0, 1.1, 0.050), minor=True)
iaxes.set_yticks(np.arange(0.0, 1.1, 0.010))
iaxes.set_yticks(np.arange(0.0, 1.1, 0.005), minor=True)
iaxes.set_xlim(0.50, 1.00)
iaxes.set_ylim(0.00, 0.05)
iaxes.grid()
#iaxes.set_aspect('equal')
ax = fig.add_subplot(122)
ax.plot(thresholds, tpr, color='tab:red')
ax.set_xlabel("Threshold", fontsize=15)
ax.set_ylabel("True Positive Rate", fontsize=15)
ax.set_ylim(-0.05,1.05)
ax.set_xlim(-0.05,1.05)
ax.tick_params(labelsize=12)
ax.set_aspect('equal')
iaxes = inset_axes(ax, width="50%", height="50%", loc=4)
ip = InsetPosition(ax, [0.15,0.10,0.7,0.3])
iaxes.set_axes_locator(ip)
iaxes.tick_params(labelsize=12)
iaxes.plot(thresholds, tpr, color='tab:red')
iaxes.set_xticks(np.arange(0.0, 1.1, 0.100))
iaxes.set_xticks(np.arange(0.0, 1.1, 0.050), minor=True)
iaxes.set_yticks(np.arange(0.0, 1.1, 0.040))
iaxes.set_yticks(np.arange(0.0, 1.1, 0.010), minor=True)
iaxes.set_xlim(0.50, 1.00)
iaxes.set_ylim(0.76, 0.96)
iaxes.grid()
#iaxes.set_aspect('equal')
#plt.savefig('TPR_FPR_vs_Threshold.png',dpi=100)
plt.show()
plt.close()
In [51]:
index = fp.index.values
print index
show_images(index[:10],*np.split(bogus_images, 3, axis=3), #labelAsFloat=True,
labels=bogus_params.index.values,
classes=bogus_params['ana_cnn_prob'],
#save_name='psstar_and_dist_less_than_3arcsec',
title='False Positive')
In [58]:
real_params['sub_snr'] = real_params['sub_flux_auto'] / real_params['sub_fluxerr_auto']
bogus_params['sub_snr'] = bogus_params['sub_flux_auto'] / bogus_params['sub_fluxerr_auto']
In [61]:
label = np.concatenate([np.ones(real_params[~((real_params['ps_type']=='star')&(real_params['ps_match']==True))&
(real_params['ana_param_cand']==True)]['ana_cnn_prob'].shape[0]),
np.zeros(bogus_params[~((bogus_params['ps_type']=='star')&(bogus_params['ps_match']==True))&
(bogus_params['ana_param_cand']==True)]['ana_cnn_prob'].shape[0])],
axis=0)
prob = np.concatenate([real_params[~((real_params['ps_type']=='star')&(real_params['ps_match']==True))&
(real_params['ana_param_cand']==True)]['ana_cnn_prob'].values,
bogus_params[~((bogus_params['ps_type']=='star')&(bogus_params['ps_match']==True))&
(bogus_params['ana_param_cand']==True)]['ana_cnn_prob'].values],
axis=0)
snr = np.concatenate([real_params[~((real_params['ps_type']=='star')&(real_params['ps_match']==True))&
(real_params['ana_param_cand']==True)]['sub_snr'].values,
bogus_params[~((bogus_params['ps_type']=='star')&(bogus_params['ps_match']==True))&
(bogus_params['ana_param_cand']==True)]['sub_snr'].values],
axis=0)
In [71]:
fig = plt.figure(figsize=(12,6))
fig.suptitle('Frequency Distribution of SNR', fontsize=18)
ax = fig.add_subplot(121)
ax.hist([bogus_params[~((bogus_params['ps_type']=='star')&(bogus_params['ps_match']==True))&(bogus_params['ana_param_cand']==True)]['sub_snr']],
label=['BOGUS Transient'],
color=['tab:red'],
range=(0,50),
bins=10, density=True)
#ax.set_xlim(0,50)
#ax.set_ylim(0,200)
ax.set_xlabel('SNR', fontsize=15)
#ax.set_yscale('log')
ax.set_ylabel('Frequency', fontsize=15)
ax.tick_params(labelsize=12)
plt.legend(fontsize=14)
ax = fig.add_subplot(122)
ax.hist([bogus_params[~((bogus_params['ps_type']=='star')&(bogus_params['ps_match']==True))&(bogus_params['ana_param_cand']==True)]['sub_snr']],
label=['BOGUS Transient'],
color=['tab:red'],
range=(0,50),
bins=10, density=True)
ax.hist([real_params[~((real_params['ps_type']=='star')&(real_params['ps_match']==True))&(real_params['ana_param_cand']==True)]['sub_snr']],
label=['REAL Transient'],
color=['tab:blue'],
range=(0,50),
bins=10, density=True)
#ax.set_xlim(0,50)
#ax.set_ylim(0,200)
ax.set_xlabel('SNR', fontsize=15)
#ax.set_yscale('log')
#ax.set_ylabel('Frequency')
ax.tick_params(labelsize=12)
plt.legend(fontsize=14)
#plt.savefig('/home/hamasaki/master_reseach/dist_snr.png')
plt.show()
plt.close()
In [72]:
def snr_roc(sn_thresh=20):
high_sn_index = np.where(snr >= sn_thresh)[0]
low_sn_index = np.where(snr < sn_thresh)[0]
fpr_high, tpr_high, thresholds_high = roc_curve(label[high_sn_index], prob[high_sn_index])
auc_high = auc(fpr_high, tpr_high)
fpr_low, tpr_low, thresholds_low = roc_curve(label[low_sn_index], prob[low_sn_index])
auc_low = auc(fpr_low, tpr_low)
return auc_high, auc_low
sn_thresh_list = range(8, 101, 1)
auc_high_list = []
auc_low_list = []
for sn_thresh in sn_thresh_list:
auc_high, auc_low = snr_roc(sn_thresh)
auc_high_list.append(auc_high)
auc_low_list.append(auc_low)
#auc_ratio = (np.array(auc_high_list) - base_auc_value)/(np.array(auc_low_list) - base_auc_value)
#print auc_high_list
#print auc_low_list
#print auc_ratio
fig = plt.figure(figsize=(15,7))
fig.suptitle('SNR vs. AUC',fontsize=18)
ax = fig.add_subplot(121)
ax.plot(np.array(sn_thresh_list), np.array(auc_high_list), '-', label='AUC (SNR > SNR threshold)', color='tab:orange')
ax.plot(np.array(sn_thresh_list), np.array(auc_low_list), '-', label='AUC (SNR < SNR threshold)', color='tab:cyan')
#ax.hlines([base_auc],8,101,'gray',linestyles='dashed')
#ax.set_xlim(0, 100)
#ax.set_ylim(0, 1500)
ax.tick_params(labelsize=12)
ax.set_xlabel('SNR threshold',fontsize=15)
ax.set_ylabel('AUC',fontsize=15)
plt.legend(fontsize=15)
#plt.show()
#plt.close()
def snr_true_num(sn_thresh=20):
high_sn_index = np.where( (snr >= sn_thresh)&(label==1) )[0]
low_sn_index = np.where( (snr < sn_thresh)&(label==1) )[0]
return [high_sn_index.shape[0], low_sn_index.shape[0]]
sn_high_low_num = [ snr_true_num(sn_thresh=thresh) for thresh in sn_thresh_list]
sn_high_num = np.array([num[0] for num in sn_high_low_num])
sn_low_num = np.array([num[1] for num in sn_high_low_num])
true_num_ratio = sn_high_num.astype('float64') / (sn_low_num + sn_high_num)
def snr_false_num(sn_thresh=20):
high_sn_index = np.where( (snr >= sn_thresh)&(label==0) )[0]
low_sn_index = np.where( (snr < sn_thresh)&(label==0) )[0]
return [high_sn_index.shape[0], low_sn_index.shape[0]]
sn_high_low_num = [ snr_false_num(sn_thresh=thresh) for thresh in sn_thresh_list]
sn_high_num = np.array([num[0] for num in sn_high_low_num])
sn_low_num = np.array([num[1] for num in sn_high_low_num])
false_num_ratio = sn_high_num.astype('float64') / (sn_low_num + sn_high_num)
#fig = plt.figure(figsize=(8,6))
#fig.suptitle('Rate of Number of High SNR VS SNR (All=1)', fontsize=20)
ax = fig.add_subplot(122)
ax.plot(np.array(sn_thresh_list), true_num_ratio, '-', label='REAL Transient', color='tab:blue')
ax.plot(np.array(sn_thresh_list), false_num_ratio, '-', label='BOGUS Transient', color='tab:red')
#ax.set_xlim(0, 100)
#ax.set_ylim(0, 1)
ax.tick_params(labelsize=12)
plt.legend(fontsize=15)
ax.set_xlabel('SNR threshold', fontsize=15)
ax.set_ylabel('Ratio of Number (SNR > SNR threshold) ', fontsize=15)
#plt.savefig('/home/hamasaki/master_reseach/snr_vs_auc.png')
plt.show()
plt.close()
In [73]:
sn_thresh = 20
high_thresh_index = np.where(snr >= sn_thresh)[0]
low_thresh_index = np.where(snr < sn_thresh)[0]
#prob_test = model0.predict_proba(image_test, verbose=False)
fpr_high, tpr_high, thresholds_high = roc_curve(label[high_thresh_index], prob[high_thresh_index])
val_auc_high = auc(fpr_high, tpr_high)
fpr_low, tpr_low, thresholds_low = roc_curve(label[low_thresh_index], prob[low_thresh_index])
val_auc_low = auc(fpr_low, tpr_low)
fig = plt.figure(figsize=(7,7))
fig.suptitle("ROC curve (AUC = %6.4f / %6.4f)" % (val_auc_high, val_auc_low), fontsize=18)
ax = fig.add_subplot(111)
ax.plot(fpr_low, tpr_low, 'k:', label='SNR < %3.1f' % (sn_thresh))
ax.plot(fpr_high, tpr_high, 'r-', label='SNR > %3.1f' % (sn_thresh))
ax.set_xlabel("False Positive Rate", fontsize=15)
ax.set_ylabel("True Positive Rate", fontsize=15)
ax.set_ylim(-0.05,1.05)
ax.set_xlim(-0.05,1.05)
#ax.tick_params(labelsize=15)
#ax.grid()
ax.set_aspect('equal')
plt.legend(bbox_to_anchor=(0.95,0.9), fontsize=14)
iaxes = inset_axes(ax, width="50%", height="50%", loc=4)
ip = InsetPosition(ax, [0.4,0.1,0.5,0.5])
iaxes.set_axes_locator(ip)
iaxes.plot(fpr_low, tpr_low, 'k:')
iaxes.plot(fpr_high, tpr_high, 'r-')
iaxes.set_xticks(np.arange(0.0, 0.2, 0.050))
iaxes.set_xticks(np.arange(0.0, 0.2, 0.010), minor=True)
iaxes.set_yticks(np.arange(0.80, 1.0, 0.020))
iaxes.set_yticks(np.arange(0.80, 1.0, 0.010), minor=True)
iaxes.set_xlim(0.00, 0.15)
iaxes.set_ylim(0.80, 0.98)
#iaxes.tick_params(labelsize=14)
iaxes.grid()
iaxes.set_aspect('equal')
plt.show()
plt.close()
