Index: /branches/alston_branches/nn.20200129/Template.py =================================================================== --- /branches/alston_branches/nn.20200129/Template.py (revision 41336) +++ /branches/alston_branches/nn.20200129/Template.py (revision 41336) @@ -0,0 +1,262 @@ +#!/usr/bin/env python2 +# -*- coding: utf-8 -*- +""" +Created on Tue Apr 7 17:23:17 2020 + +@author: jack +""" + +import numpy as np +import matplotlib.pyplot as plt +import Pan_STARRS_galaxies_final as ps +from astropy.io import fits +import network + +#specify filepaths + +smf_filepath = +raw_data_loacation = +destination_location = +print('smf data specified') + +#plot kron-inst against inst to check cuts + +def data_extractor_all(filepath): + hdu_list = fits.open(filepath) + data = [hdu_list[3*i+3].data for i in range(len(hdu_list)/3)] + names = [hdu_list[3*i+3].header[-1] for i in range(len(hdu_list)/3)] + empty_array = np.ndarray.tolist(np.zeros(len(data))) + for i in range(len(data)): + element = np.ndarray.tolist(np.zeros(4)) + element[0] = names[i] + element[1] = [hdu_list[3*i+3].data[j][0] for j in range(len(hdu_list[3*i+3].data))] #object id + element[2] = [hdu_list[3*i+3].data[j][7] for j in range(len(hdu_list[3*i+3].data))] #inst mag + element[3] = [-2.5*np.log10(hdu_list[3*i+3].data[j][46]) for j in range(len(hdu_list[3*i+3].data))] #kron flux + empty_array[i] = element + return empty_array + +data = data_extractor_all(smf_filepath) +kron = [] +inst = [] +for i in range(len(data)): + for j in range(len(data[i][3])): + kron.append(data[i][3][j]) + inst.append(data[i][2][j]) +diff = np.array(kron)- np.array(inst) +fig = plt.figure(figsize = [10,10]) +ax = plt.gca() +ax.plot(inst, diff, ',') +ax.set_xlim(-17.5, -9.5) +ax.set_ylim(-2, 2) + +mag_cut = -9.5 + +#extract stellar objects and masks + +def data_extractor_MS(filepath): + hdu_list = fits.open(filepath) + data = [hdu_list[3*i+3].data for i in range(len(hdu_list)/3)] + names = [hdu_list[3*i+3].header[-1] for i in range(len(hdu_list)/3)] + empty_array = np.ndarray.tolist(np.zeros(len(data))) + for i in range(len(data)): + element = np.ndarray.tolist(np.zeros(8)) + psf_qf = [hdu_list[3*i+3].data[j][33] for j in range(len(hdu_list[3*i+3].data))] + inst_mag = [hdu_list[3*i+3].data[j][7] for j in range(len(hdu_list[3*i+3].data))] + difference = [-2.5*np.log10(hdu_list[3*i+3].data[j][46])-hdu_list[3*i+3].data[j][7] for j in range(len(hdu_list[3*i+3].data))] + kron_flux = [hdu_list[3*i+3].data[j][46] for j in range(len(hdu_list[3*i+3].data))] + print(i) + element[0] = names[i] + element[1] = [hdu_list[3*i+3].data[j][0] for j in range(len(hdu_list[3*i+3].data)) if psf_qf[j] > 0.995 and str(kron_flux[j]) != 'nan' and difference[j] > -0.2 and inst_mag[j] < mag_cut] #object id + element[2] = [hdu_list[3*i+3].data[j][7] for j in range(len(hdu_list[3*i+3].data)) if psf_qf[j] > 0.995 and str(kron_flux[j]) != 'nan' and difference[j] > -0.2 and inst_mag[j] < mag_cut] #inst mag + element[3] = [hdu_list[3*i+3].data[j][8] for j in range(len(hdu_list[3*i+3].data)) if psf_qf[j] > 0.995 and str(kron_flux[j]) != 'nan' and difference[j] > -0.2 and inst_mag[j] < mag_cut] #inst mag sig + element[4] = [hdu_list[3*i+3].data[j][33] for j in range(len(hdu_list[3*i+3].data)) if psf_qf[j] > 0.995 and str(kron_flux[j]) != 'nan' and difference[j] > -0.2 and inst_mag[j] < mag_cut] #psf qf + element[5] = [hdu_list[3*i+3].data[j][34] for j in range(len(hdu_list[3*i+3].data)) if psf_qf[j] > 0.995 and str(kron_flux[j]) != 'nan' and difference[j] > -0.2 and inst_mag[j] < mag_cut] #psf qf perfect + element[6] = [-2.5*np.log10(hdu_list[3*i+3].data[j][46]) for j in range(len(hdu_list[3*i+3].data)) if psf_qf[j] > 0.995 and str(kron_flux[j]) != 'nan' and difference[j] > -0.2 and inst_mag[j] < mag_cut] #kron flux + element[7] = [hdu_list[3*i+3].data[j][53] & 0x2000 for j in range(len(hdu_list[3*i+3].data)) if psf_qf[j] > 0.995 and str(kron_flux[j]) != 'nan' and difference[j] > -0.2 and inst_mag[j] < mag_cut] #checks if mask is present + empty_array[i] = element + return empty_array + +SM_data = data_extractor_MS(smf_filepath) +print('SM data extracted') + +#save stellar objects and masks + +ps.file_saver_MS(SM_data, raw_data_loacation, destination_location) +print('SM data saved') + +#extract extended objects + +def data_extractor_E(filepath): + hdu_list = fits.open(filepath) + data = [hdu_list[3*i+3].data for i in range(len(hdu_list)/3)] + names = [hdu_list[3*i+3].header[-1] for i in range(len(hdu_list)/3)] + empty_array = np.ndarray.tolist(np.zeros(len(data))) + for i in range(len(data)): + element = np.ndarray.tolist(np.zeros(8)) + psf = [hdu_list[3*i+3].data[j][33] for j in range(len(hdu_list[3*i+3].data))] + mask_present = [hdu_list[3*i+3].data[j][53] & 0x2000 for j in range(len(hdu_list[3*i+3].data))] #whack hexadecimal shit + saturated = [hdu_list[3*i+3].data[j][53] & 0x80 for j in range(len(hdu_list[3*i+3].data))] + kron_flux = [hdu_list[3*i+3].data[j][46] for j in range(len(hdu_list[3*i+3].data))] + difference = [-2.5*np.log10(hdu_list[3*i+3].data[j][46])-hdu_list[3*i+3].data[j][7] for j in range(len(hdu_list[3*i+3].data))] + mag = [hdu_list[3*i+3].data[j][7] for j in range(len(hdu_list[3*i+3].data))] + print(i) + element[0] = names[i] + element[1] = [hdu_list[3*i+3].data[j][0] for j in range(len(hdu_list[3*i+3].data)) if mask_present[j] == 0 and str(kron_flux[j]) != 'nan' and psf[j] > 0.995 and difference[j] < -0.2 and saturated[j] == 0 and mag[j] < mag_cut] #object id + element[2] = [hdu_list[3*i+3].data[j][7] for j in range(len(hdu_list[3*i+3].data)) if mask_present[j] == 0 and str(kron_flux[j]) != 'nan' and psf[j] > 0.995 and difference[j] < -0.2 and saturated[j] == 0 and mag[j] < mag_cut]#inst mag + element[3] = [hdu_list[3*i+3].data[j][8] for j in range(len(hdu_list[3*i+3].data)) if mask_present[j] == 0 and str(kron_flux[j]) != 'nan' and psf[j] > 0.995 and difference[j] < -0.2 and saturated[j] == 0 and mag[j] < mag_cut] #inst mag sig + element[4] = [hdu_list[3*i+3].data[j][33] for j in range(len(hdu_list[3*i+3].data)) if mask_present[j] == 0 and str(kron_flux[j]) != 'nan' and psf[j] > 0.995 and difference[j] < -0.2 and saturated[j] == 0 and mag[j] < mag_cut] #psf qf + element[5] = [hdu_list[3*i+3].data[j][34] for j in range(len(hdu_list[3*i+3].data)) if mask_present[j] == 0 and str(kron_flux[j]) != 'nan' and psf[j] > 0.995 and difference[j] < -0.2 and saturated[j] == 0 and mag[j] < mag_cut] #psf qf perfect + element[6] = [-2.5*np.log10(hdu_list[3*i+3].data[j][46]) for j in range(len(hdu_list[3*i+3].data)) if mask_present[j] == 0 and str(kron_flux[j]) != 'nan' and psf[j] > 0.995 and difference[j] < -0.2 and saturated[i] == 0 and mag[j] < mag_cut] #kron flux + element[7] = [hdu_list[3*i+3].data[j][53] for j in range(len(hdu_list[3*i+3].data)) if mask_present[j] == 0 and str(kron_flux[j]) != 'nan' and psf[j] > 0.995 and mag[j] < -9.5 and difference[j] < -0.2 and saturated[j] == 0 and mag[j] < mag_cut] #flags + empty_array[i] = element + return empty_array + +E_data = data_extractor_E(smf_filepath) +print('E data extracted') + +#manually sort extended objects + +ps.file_saver_E(E_data, raw_data_loacation, destination_location, 5) +print('E data saved') + +#generate testing data + +trn_filenames_S = ps.filenames(destination_location + r'/training/star/') +trn_filenames_E = ps.filenames(destination_location + r'/training/extended good/') +trn_filenames_M = ps.filenames(destination_location + r'/training/masked/') + +tst_filenames_S = ps.filenames(destination_location + r'/testing/star/') +tst_filenames_E = ps.filenames(destination_location + r'/testing/extended good/') +tst_filenames_M = ps.filenames(destination_location + r'/testing/masked/') + +testing_filenames = trn_filenames_S + tst_filenames_S + trn_filenames_E + tst_filenames_E + trn_filenames_M + tst_filenames_M + +testing_MO = ps.testing_set([trn_filenames_S + tst_filenames_S + trn_filenames_E + tst_filenames_E, trn_filenames_M + tst_filenames_M], [1,1]) +testing_SE = ps.testing_set([trn_filenames_S + tst_filenames_S, trn_filenames_E + tst_filenames_E], [1,1]) + +'''testing_MO = np.load(, allow_pickle = True) +testing_SE = np.load(, allow_pickle = True)''' #may load data using this if testing sets already generated + +print('testing data generated') + +#feed data through neural networks + +sizes_SE = np.load(r'/home/jack/Documents/Research Project/Main code/star galaxy mask/Network Data/magnitude limited networks/deep learning test/best_array_sizes.npy', allow_pickle = True) +weights_SE = np.load(r'/home/jack/Documents/Research Project/Main code/star galaxy mask/Network Data/magnitude limited networks/deep learning test/best_array_weights.npy', allow_pickle = True) +biases_SE = np.load(r'/home/jack/Documents/Research Project/Main code/star galaxy mask/Network Data/magnitude limited networks/deep learning test/best_array_biases.npy', allow_pickle = True) +net_SE = network.Network(sizes_SE) +net_SE.weights = weights_SE +net_SE.biases = biases_SE + +sizes_MO = np.load(r'/home/jack/Documents/Research Project/Main code/star galaxy mask/Network Data/magnitude limited networks/mask other/best_array_sizes.npy', allow_pickle = True) +weights_MO = np.load(r'/home/jack/Documents/Research Project/Main code/star galaxy mask/Network Data/magnitude limited networks/mask other/best_array_weights.npy', allow_pickle = True) +biases_MO = np.load(r'/home/jack/Documents/Research Project/Main code/star galaxy mask/Network Data/magnitude limited networks/mask other/best_array_biases.npy', allow_pickle = True) +net_MO = network.Network(sizes_MO) +net_MO.weights = weights_MO +net_MO.biases = biases_MO +print('networks loaded') + +#plot neural network results +#note that string indices eg. [99:107] will need to be changed to suit varying file paths. + +useful_results = ps.network_results(smf_filepath, testing_filenames) +'''useful_results = np.load(, allow_pickle = True)''' #if results already extracted, loads file + +chip_id = [useful_results[0][i] for i in range(len(useful_results[0])) if float(useful_results[3][i]) < -9.5] +image_id = [useful_results[1][i] for i in range(len(useful_results[0])) if float(useful_results[3][i]) < -9.5] +kron_mag = [float(useful_results[2][i]) for i in range(len(useful_results[0])) if float(useful_results[3][i]) < -9.5] +inst_mag = [float(useful_results[3][i]) for i in range(len(useful_results[0])) if float(useful_results[3][i]) < -9.5] +diff = [float(useful_results[4][i]) for i in range(len(useful_results[0])) if float(useful_results[3][i]) < -9.5] +testing_fnames = [testing_filenames[i] for i in range(len(useful_results[0])) if float(useful_results[3][i]) < -9.5] + +#run data through networks + +mask = [] +other = [] +for i in range(len(testing_MO)): + result = np.argmax(net_MO.feedforward(testing_MO[i][0])) + network_output = net_MO.feedforward(testing_MO[i][0]) + if result == 0: + other.append([testing_MO[i][0], chip_id[i], image_id[i], kron_mag[i], inst_mag[i], diff[i], testing_fnames[i], network_output]) + if result == 1: + mask.append([testing_MO[i][0], chip_id[i], image_id[i], kron_mag[i], inst_mag[i], diff[i], testing_fnames[i], network_output]) + +star = [] +extended = [] +for i in range(len(other)): + result = np.argmax(net_SE.feedforward(other[i][0])) + network_output = net_SE.feedforward(other[i][0]) + if result == 0: + star.append([other[i][0], other[i][1], other[i][2], other[i][3], other[i][4], other[i][5], other[i][6], network_output]) + if result == 1: + extended.append([other[i][0], other[i][1], other[i][2], other[i][3], other[i][4], other[i][5], other[i][6], network_output]) + +#Anyalyse results +no_extended_right = float(len([i[6] for i in extended if i[6][99:107] == 'extended'])) +no_star_right = float(len([i[6] for i in star if i[6][99:103] == 'star'])) +no_mask_right = float(len([i[6] for i in mask if i[6][99:105] == 'masked'])) +no_extended = float(len([i for i in testing_SE if i[1] == 1])) +no_star = float(len([i for i in testing_SE if i[1] == 0])) +no_mask = float(len([i for i in testing_MO if i[1] == 1])) +print(100*no_star_right/no_star, 100*no_mask_right/no_mask, 100*no_extended_right/no_extended) #print percentage effectiveness + + +star_right_inst = [i[4] for i in star if i[6][99:103] == 'star'] +star_wrong_extended_inst = [i[4] for i in star if i[6][99:107] == 'extended'] +star_wrong_mask_inst = [i[4] for i in star if i[6][99:105] == 'masked'] + +extended_right_inst = [i[4] for i in extended if i[6][99:107] == 'extended'] +extended_wrong_mask_inst = [i[4] for i in extended if i[6][99:105] == 'masked'] +extended_wrong_star_inst = [i[4] for i in extended if i[6][99:103] == 'star'] + +mask_right_inst = [i[4] for i in mask if i[6][99:105] == 'masked'] +mask_wrong_star_inst = [i[4] for i in mask if i[6][99:103] == 'star'] +mask_wrong_extended_inst = [i[4] for i in mask if i[6][99:107] == 'extended'] + +star_right_kron = [i[5] for i in star if i[6][99:103] == 'star'] +star_wrong_extended_kron = [i[5] for i in star if i[6][99:107] == 'extended'] +star_wrong_mask_kron = [i[5] for i in star if i[6][99:105] == 'masked'] + +extended_right_kron = [i[5] for i in extended if i[6][99:107] == 'extended'] +extended_wrong_mask_kron = [i[5] for i in extended if i[6][99:105] == 'masked'] +extended_wrong_star_kron = [i[5] for i in extended if i[6][99:103] == 'star'] + +mask_right_kron = [i[5] for i in mask if i[6][99:105] == 'masked'] +mask_wrong_star_kron = [i[5] for i in mask if i[6][99:103] == 'star'] +mask_wrong_extended_kron = [i[5] for i in mask if i[6][99:107] == 'extended'] + +fig1 = plt.figure(figsize=[10,10]) +ax1 = fig1.gca() + +ax1.plot(star_right_inst, star_right_kron, 'r.', label = 'Correctly Identified Point Source') +ax1.plot(extended_wrong_star_inst, extended_wrong_star_kron, 'g.', label = 'Point Source Identified as Extended Source') +ax1.plot(mask_wrong_star_inst, mask_wrong_star_kron, 'b.', label = 'Point Source Identified as Masked Ray') + +ax1.plot(star_wrong_mask_inst, star_wrong_mask_kron, 'rx', label = 'Masked Ray Identified as Point Source') +ax1.plot(mask_right_inst, mask_right_kron, 'bx', label = 'Correctly Identified Masked Ray') +ax1.plot(extended_wrong_mask_inst, extended_wrong_mask_kron, 'gx', label = 'Masked Ray Identified as Extended Source') + +ax1.scatter(star_wrong_extended_inst, star_wrong_extended_kron, edgecolors = 'r', facecolors = 'none', label = ' Extended source Identified as Point Source') +ax1.scatter(extended_right_inst, extended_right_kron, edgecolors = 'g', facecolors = 'none', label = 'Correctly Identified Extended Source') +ax1.scatter(mask_wrong_extended_inst, mask_wrong_extended_kron, edgecolors = 'b', facecolors = 'none', label = 'Extended Source Identified as Masked Ray') + +ax1.legend() +ax1.set_xlabel('Instrumental Magnitude') +ax1.set_ylabel('Kron Magnitude - Intrusmental Magnitude') + + +star_right_fnames = [i[6] for i in star if i[6][81:85] == 'star'] +star_wrong_extended_fnames = [i[6] for i in star if i[6][81:89] == 'extended'] +star_wrong_mask_fnames = [i[6] for i in star if i[6][81:87] == 'masked'] + +extended_right_fnames = [i[6] for i in extended if i[6][81:89] == 'extended'] +extended_wrong_mask_fnames = [i[6] for i in extended if i[6][81:87] == 'masked'] +extended_wrong_star_fnames = [i[6] for i in extended if i[6][81:85] == 'star'] + +mask_right_fnames = [i[6] for i in mask if i[6][81:87] == 'masked'] +mask_wrong_star_fnames = [i[6] for i in mask if i[6][81:85] == 'star'] +mask_wrong_extended_fnames = [i[6] for i in mask if i[6][81:89] == 'extended'] + +print(len(star_right_fnames), len(star_wrong_extended_fnames), len(star_wrong_mask_fnames)) +print(len(extended_wrong_star_fnames), len(extended_right_fnames), len(extended_wrong_mask_fnames)) +print(len(mask_wrong_star_fnames), len(mask_wrong_extended_fnames), len(mask_right_fnames))