nustar/nuwavdet
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added Det2raw for region files, PEP8 adjustments

Browse Source
This commit is contained in:
Andrey Mukhin 2022-09-19 12:31:33 +03:00
parent 8db10f4bec
commit 094063d21d
1 changed files with 212 additions and 133 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

345
nuwavsource/nuwavsource.py
View File

@ -2,7 +2,7 @@
import numpy as np import numpy as np
import itertools import itertools
from pandas import DataFrame, read_csv from pandas import DataFrame, read_csv
from astropy.table import Table from astropy.table import Table, unique
from astropy.coordinates import SkyCoord from astropy.coordinates import SkyCoord
from astropy import units as u from astropy import units as u
from multiprocessing import get_context, cpu_count from multiprocessing import get_context, cpu_count
@ -16,25 +16,38 @@ from glob import glob
from warnings import filterwarnings from warnings import filterwarnings
filterwarnings('ignore') filterwarnings('ignore')
def get_link_list(folder, sort_list=True): def get_link_list(folder, sort_list=True):
links = glob(f'{folder}\\**\\*_cl.evt',recursive=True) links = glob(f'{folder}\\**\\*_cl.evt', recursive=True)
if sort_list: if sort_list:
sorted_list = sorted(links, key=lambda x: stat(x).st_size) sorted_list = sorted(links, key=lambda x: stat(x).st_size)
return np.array(sorted_list) return np.array(sorted_list)
else: else:
return np.array(links) return np.array(links)
def binary_array(num): def binary_array(num):
variants = [[0,1],]*num variants = [[0, 1], ]*num
out = np.zeros((2**num,num),bool) out = np.zeros((2**num, num), bool)
for idx, level in enumerate(itertools.product(*variants)): for idx, level in enumerate(itertools.product(*variants)):
out[idx] = np.array(level,dtype=bool) out[idx] = np.array(level, dtype=bool)
return out return out
def create_array(filename,mode='Sky'):
temp = fits.getdata(filename,1)
def create_array(filename, mode='Sky'):
temp = fits.getdata(filename, 1)
if mode == 'Sky': if mode == 'Sky':
return np.histogram2d(temp['Y'],temp['X'],1000,[[0, 1000], [0, 1000]])[0] return np.histogram2d(temp['Y'],
temp['X'],
1000,
[[0, 1000], [0, 1000]])[0]
if mode == 'Det': if mode == 'Det':
return np.histogram2d(temp['DET1Y'],temp['DET1X'],360,[[0, 360], [0, 360]])[0] return np.histogram2d(temp['DET1Y'],
temp['DET1X'],
360,
[[0, 360], [0, 360]])[0]
def get_wcs(file): def get_wcs(file):
header = file[1].header header = file[1].header
wcs = WCS({ wcs = WCS({
@ -46,7 +59,9 @@ def get_wcs(file):
'NAXIS1': header['TLMAX38'], 'NAXIS2': header['TLMAX39'] 'NAXIS1': header['TLMAX38'], 'NAXIS2': header['TLMAX39']
}) })
return wcs return wcs
def atrous(level = 0, max_size = 1001):
def atrous(level=0, max_size=1001):
base = 1/256*np.array([ base = 1/256*np.array([
[1, 4, 6, 4, 1], [1, 4, 6, 4, 1],
[4, 16, 24, 16, 4], [4, 16, 24, 16, 4],
@ -55,67 +70,82 @@ def atrous(level = 0, max_size = 1001):
[1, 4, 6, 4, 1], [1, 4, 6, 4, 1],
]) ])
size = 2**level * (base.shape[0]-1)+1 size = 2**level * (base.shape[0]-1)+1
output = np.zeros((size,size)) output = np.zeros((size, size))
output[::2**level, ::2**level] = base output[::2**level, ::2**level] = base
if output.shape[0]>max_size: if output.shape[0] > max_size:
return output[(size-1)//2-(max_size-1)//2:(size-1)//2+(max_size-1)//2+1, (size-1)//2-(max_size-1)//2:(size-1)//2+(max_size-1)//2+1] return output[(size-1)//2-(max_size-1)//2:(size-1)//2+(max_size-1)//2+1,
(size-1)//2-(max_size-1)//2:(size-1)//2+(max_size-1)//2+1]
return output return output
def gauss(level=0, max_size = 1000):
def gauss(level=0, max_size=1000):
size = min(5*2**(level+1)+1, max_size) size = min(5*2**(level+1)+1, max_size)
sigma = 2**(level) sigma = 2**(level)
A = 1/(2*np.pi*sigma**2)**0.5 A = 1/(2*np.pi*sigma**2)**0.5
x = A*np.exp( (-(np.arange(size)-(size-1)//2)**2)/(2*sigma**2)) x = A*np.exp((-(np.arange(size)-(size-1)//2)**2)/(2*sigma**2))
out = np.multiply.outer(x,x) out = np.multiply.outer(x, x)
return out return out
def adjecent(array): def adjecent(array):
grid = np.array([ grid = np.array([
[1,1,1], [1, 1, 1],
[1,0,1], [1, 0, 1],
[1,1,1] [1, 1, 1]
]) ])
output = fftconvolve(array,grid,mode='same') >= 0.5 output = fftconvolve(array, grid, mode='same') >= 0.5
try: try:
output = np.logical_and(np.logical_and(output, np.logical_not(array)),np.logical_not(array.mask)) output = np.logical_and(np.logical_and(output, np.logical_not(array)),
np.logical_not(array.mask))
except AttributeError: except AttributeError:
output = np.logical_and(output, np.logical_not(array)) output = np.logical_and(output, np.logical_not(array))
output = np.argwhere(output == True) output = np.argwhere(output == True)
return output[:,0], output[:,1] return output[:, 0], output[:, 1]
def add_borders(array,middle=True):
mask = np.zeros(array.shape)
datax, datay = np.any(array>0,0), np.any(array>0,1) def add_borders(array, middle=True):
#Add border masks mask = np.zeros(array.shape)
x_min, y_min = np.argmax(datax), np.argmax(datay) datax, datay = np.any(array > 0, 0), np.any(array > 0, 1)
x_max, y_max = len(datax) - np.argmax(datax[::-1]), len(datay) - np.argmax(datay[::-1]) # Add border masks
mask[y_min:y_max,x_min:x_max] = True x_min, y_min = np.argmax(datax), np.argmax(datay)
if middle is True: x_max = len(datax) - np.argmax(datax[::-1])
mask[176:191,:] = False y_max = len(datay) - np.argmax(datay[::-1])
mask[:,176:191] = False mask[y_min:y_max, x_min:x_max] = True
mask = np.logical_not(mask) if middle is True:
return mask mask[176:191, :] = False
mask[:, 176:191] = False
mask = np.logical_not(mask)
return mask
def fill_poisson(array, size_input=32): def fill_poisson(array, size_input=32):
if not(isinstance(array,np.ma.MaskedArray)): if not (isinstance(array, np.ma.MaskedArray)):
print('No mask found') print('No mask found')
return array return array
size = size_input size = size_input
output = array.data.copy() output = array.data.copy()
mask = array.mask.copy() mask = array.mask.copy()
while mask.sum()>1: while mask.sum() > 1:
kernel = np.ones((size,size))/size**2 kernel = np.ones((size, size))/size**2
coeff = fftconvolve(np.logical_not(mask),kernel,mode='same') coeff = fftconvolve(np.logical_not(mask), kernel, mode='same')
mean = fftconvolve(output,kernel,mode='same') mean = fftconvolve(output, kernel, mode='same')
idx = np.where(np.logical_and(mask,coeff>0.1)) idx = np.where(np.logical_and(mask, coeff > 0.1))
output[idx] = np.random.poisson(np.abs(mean[idx]/coeff[idx])) output[idx] = np.random.poisson(np.abs(mean[idx]/coeff[idx]))
mask[idx] = False mask[idx] = False
size *= 2 size *= 2
return output return output
def mirror(array): def mirror(array):
size = array.shape[0] size = array.shape[0]
output = np.tile(array,(3,3)) output = np.tile(array, (3, 3))
output[0:size] = np.flipud(output[0:size]) output[0:size] = np.flipud(output[0:size])
output[2*size:3*size] = np.flipud(output[2*size:3*size]) output[2*size:3*size] = np.flipud(output[2*size:3*size])
output[:,0:size] = np.fliplr(output[:,0:size]) output[:, 0:size] = np.fliplr(output[:, 0:size])
output[:,2*size:3*size] = np.fliplr(output[:,2*size:3*size]) output[:, 2*size:3*size] = np.fliplr(output[:, 2*size:3*size])
return output return output
class Observation: class Observation:
def __init__(self, file_name, E_borders=[3,20]): def __init__(self, file_name, E_borders=[3,20]):
self.filename = file_name self.filename = file_name
@ -128,120 +158,154 @@ class Observation:
self.header = file[0].header self.header = file[0].header
self.det = self.header['INSTRUME'][-1] self.det = self.header['INSTRUME'][-1]
self.wcs = get_wcs(file) self.wcs = get_wcs(file)
self.data = np.ma.masked_array(*self.get_data(file,E_borders)) self.data = np.ma.masked_array(*self.get_data(file, E_borders))
self.hard_mask = add_borders(self.data.data,middle=False) self.hard_mask = add_borders(self.data.data, middle=False)
self.exposure = self.header['EXPOSURE'] self.exposure = self.header['EXPOSURE']
def get_coeff(self): def get_coeff(self):
coeff = np.array([0.977,0.861,1.163,1.05]) if self.det=='A' else np.array([1.004,0.997,1.025,0.979]) coeff = np.array([0.977, 0.861, 1.163, 1.05]) if self.det == 'A' else np.array([1.004, 0.997, 1.025, 0.979])
resized_coeff = (coeff).reshape(2,2).repeat(180,0).repeat(180,1) resized_coeff = (coeff).reshape(2, 2).repeat(180, 0).repeat(180, 1)
return resized_coeff return resized_coeff
def get_data(self, file, E_borders=[3,20]):
def get_data(self, file, E_borders=[3, 20]):
PI_min, PI_max = (np.array(E_borders)-1.6)/0.04 PI_min, PI_max = (np.array(E_borders)-1.6)/0.04
data = file[1].data.copy() data = file[1].data.copy()
idx_mask = (data['STATUS'].sum(1) == 0) idx_mask = (data['STATUS'].sum(1) == 0)
idx_output = np.logical_and(idx_mask,np.logical_and((data['PI']>PI_min),(data['PI']<PI_max))) idx_output = np.logical_and(idx_mask, np.logical_and((data['PI'] > PI_min), (data['PI'] < PI_max)))
data_output = data[idx_output] data_output = data[idx_output]
data_mask = data[np.logical_not(idx_mask)] data_mask = data[np.logical_not(idx_mask)]
build_hist = lambda array: np.histogram2d(array['DET1Y'],array['DET1X'],360,[[0,360],[0,360]])[0] build_hist = lambda array: np.histogram2d(array['DET1Y'], array['DET1X'], 360, [[0, 360], [0, 360]])[0]
output = build_hist(data_output) output = build_hist(data_output)
mask = build_hist(data_mask) mask = build_hist(data_mask)
mask = np.logical_or(mask,add_borders(output)) mask = np.logical_or(mask, add_borders(output))
mask = np.logical_or(mask, self.get_bad_pix(file)) mask = np.logical_or(mask, self.get_bad_pix(file))
return output, mask return output, mask
def get_bad_pix(self, file): def get_bad_pix(self, file):
output = np.zeros((360,360)) output = np.zeros((360, 360))
kernel = np.ones((5,5)) kernel = np.ones((5, 5))
for i in range(4): for i in range(4):
current_dir = dirname(__file__) current_dir = dirname(__file__)
pixpos_file = fits.getdata(f'{current_dir}\\pixpos\\nu{self.det}pixpos20100101v007.fits',i+1) pixpos_file = fits.getdata(f'{current_dir}\\pixpos\\nu{self.det}pixpos20100101v007.fits',i+1)
bad_pix_file = file[3+i].data.copy() bad_pix_file = file[3+i].data.copy()
temp = np.zeros(len(pixpos_file),dtype=bool) temp = np.zeros(len(pixpos_file), dtype=bool)
for x,y in zip(bad_pix_file['rawx'],bad_pix_file['rawy']): for x, y in zip(bad_pix_file['rawx'], bad_pix_file['rawy']):
temp = np.logical_or(temp, np.equal(pixpos_file['rawx'],x)*np.equal(pixpos_file['rawy'],y)) temp = np.logical_or(temp, np.equal(pixpos_file['rawx'], x)*np.equal(pixpos_file['rawy'], y))
temp = pixpos_file[temp] temp = pixpos_file[temp]
output += np.histogram2d(temp['REF_DET1Y'],temp['REF_DET1X'], 360, [[0,360],[0,360]])[0] output += np.histogram2d(temp['REF_DET1Y'], temp['REF_DET1X'], 360, [[0, 360],[0, 360]])[0]
output = convolve2d(output, kernel, mode='same') > 0 output = convolve2d(output, kernel, mode='same') > 0
return output return output
def wavdecomp(self, mode = 'gauss', thresh=False,occ_coeff = False):
#THRESHOLD def wavdecomp(self, mode='gauss', thresh=False, occ_coeff=False):
if type(thresh) is int: thresh_max, thresh_add = thresh, thresh/2 # THRESHOLD
elif type(thresh) is tuple: thresh_max, thresh_add = thresh if type(thresh) is int:
#INIT NEEDED VARIABLES thresh_max, thresh_add = thresh, thresh/2
elif type(thresh) is tuple:
thresh_max, thresh_add = thresh
# INIT NEEDED VARIABLES
wavelet = globals()[mode] wavelet = globals()[mode]
max_level = 8 max_level = 8
conv_out = np.zeros((max_level+1,self.data.shape[0],self.data.shape[1])) conv_out = np.zeros((max_level+1, self.data.shape[0], self.data.shape[1]))
size = self.data.shape[0] size = self.data.shape[0]
#PREPARE ORIGINAL DATA FOR ANALYSIS: FILL THE HOLES + MIRROR + DETECTOR CORRECTION # PREPARE ORIGINAL DATA FOR ANALYSIS: FILL THE HOLES + MIRROR + DETECTOR CORRECTION
data = fill_poisson(self.data) data = fill_poisson(self.data)
if occ_coeff: data = data*self.get_coeff() if occ_coeff:
data = data*self.get_coeff()
data = mirror(data) data = mirror(data)
data_bkg = data.copy() data_bkg = data.copy()
#ITERATIVELY CONDUCT WAVLET DECOMPOSITION # ITERATIVELY CONDUCT WAVLET DECOMPOSITION
for i in range(max_level): for i in range(max_level):
conv = fftconvolve(data,wavelet(i),mode='same') conv = fftconvolve(data, wavelet(i), mode='same')
temp_out = data-conv temp_out = data-conv
#ERRORMAP CALCULATION # ERRORMAP CALCULATION
if thresh_max != 0: if thresh_max != 0:
sig = ((wavelet(i)**2).sum())**0.5 sig = ((wavelet(i)**2).sum())**0.5
bkg = fftconvolve(data_bkg, wavelet(i),mode='same') bkg = fftconvolve(data_bkg, wavelet(i), mode='same')
bkg[bkg<0] = 0 bkg[bkg < 0] = 0
# err = (1+np.sqrt(bkg/sig**2 + 0.75))*sig**3
err = (1+np.sqrt(bkg+0.75))*sig err = (1+np.sqrt(bkg+0.75))*sig
# significant = (np.abs(temp_out)> thresh_max*err)[size:2*size,size:2*size] significant = (temp_out > thresh_max*err)[size:2*size, size:2*size]
significant = (temp_out > thresh_max*err)[size:2*size,size:2*size]
if thresh_add != 0: if thresh_add != 0:
# add_significant = (np.abs(temp_out)> thresh_add*err)[size:2*size,size:2*size] add_significant = (temp_out > thresh_add*err)[size:2*size, size:2*size]
add_significant = (temp_out > thresh_add*err)[size:2*size,size:2*size]
adj = adjecent(significant) adj = adjecent(significant)
add_condition = np.logical_and(add_significant[adj[0],adj[1]],np.logical_not(significant[adj[0],adj[1]])) add_condition = np.logical_and(add_significant[adj[0], adj[1]],
np.logical_not(significant[adj[0], adj[1]]))
while (add_condition).any(): while (add_condition).any():
to_add = adj[0][add_condition], adj[1][add_condition] to_add = adj[0][add_condition], adj[1][add_condition]
significant[to_add[0],to_add[1]] = True significant[to_add[0], to_add[1]] = True
adj = adjecent(significant) adj = adjecent(significant)
add_condition = np.logical_and(add_significant[adj[0],adj[1]],np.logical_not(significant[adj[0],adj[1]])) add_condition = np.logical_and(add_significant[adj[0], adj[1]],
# add_condition = np.logical_and(np.abs(temp_out)[adj[0],adj[1]] >= thresh_add*err[adj[0],adj[1]], np.logical_not(significant)[adj[0],adj[1]]) np.logical_not(significant[adj[0],adj[1]]))
temp_out[size:2*size,size:2*size][np.logical_not(significant)] = 0 temp_out[size:2*size, size:2*size][np.logical_not(significant)] = 0
#WRITING THE WAVELET DECOMP LAYER # WRITING THE WAVELET DECOMP LAYER
conv_out[i] = +temp_out[size:2*size,size:2*size] conv_out[i] = +temp_out[size:2*size, size:2*size]
conv_out[i][conv_out[i]<0]=0 #leave only positive data to prevent problems while summing layers # DISCARDING NEGATIVE COMPONENTS OF WAVELETS TO MAKE MASK BY SUMMING WAVELET LAYERS
conv_out[i][conv_out[i] < 0] = 0
data = conv data = conv
conv_out[max_level] = conv[size:2*size,size:2*size] conv_out[max_level] = conv[size:2*size, size:2*size]
return conv_out return conv_out
def process(obs_path):
def region_to_raw(self, region):
x_region, y_region = np.where(region)
tables = []
for i in range(4):
current_dir = dirname(__file__)
pixpos = Table(fits.getdata(f'{current_dir}\\pixpos\\nu{self.det}pixpos20100101v007.fits', i+1))
pixpos = pixpos[pixpos['REF_DET1X'] != -1]
test = np.zeros(len(pixpos['REF_DET1X']), dtype=bool)
for idx, (x, y) in enumerate(zip(pixpos['REF_DET1X'], pixpos['REF_DET1Y'])):
test[idx] = np.logical_and(np.equal(x, x_region), np.equal(y, y_region)).any()
tables.append(unique(Table({'RAWX': pixpos['RAWX'][test], 'RAWY': pixpos['RAWY'][test]})))
hdu_list = fits.HDUList([
fits.PrimaryHDU(),
fits.table_to_hdu(tables[0]),
fits.table_to_hdu(tables[1]),
fits.table_to_hdu(tables[2]),
fits.table_to_hdu(tables[3]),
])
return hdu_list
def process(args):
"""
Creates a mask using wavelet decomposition and produces some statistical and metadata about the passed observation.
args must contain two arguments: path to the file of interest and threshold, e.g. ('D:\Data\obs_cl.evt',(5,2))
"""
obs_path, thresh = args
bin_num = 6 bin_num = 6
try: try:
obs = Observation(obs_path) obs = Observation(obs_path)
sky_coord = SkyCoord(ra=obs.ra*u.deg,dec=obs.dec*u.deg,frame='fk5').transform_to('galactic') sky_coord = SkyCoord(ra=obs.ra*u.deg, dec=obs.dec*u.deg, frame='fk5').transform_to('galactic')
lon, lat = sky_coord.l.value, sky_coord.b.value lon, lat = sky_coord.l.value, sky_coord.b.value
rem_signal, rem_area, poiss_comp, rms = np.zeros((4,2**bin_num)) rem_signal, rem_area, poiss_comp, rms = np.zeros((4, 2**bin_num))
region = np.zeros(obs.data.shape, dtype=bool) region = np.zeros(obs.data.shape, dtype=bool)
rem_region = np.logical_and(region, np.logical_not(obs.data.mask)) rem_region = np.logical_and(region, np.logical_not(obs.data.mask))
masked_obs = np.ma.masked_array(obs.data, mask = region) masked_obs = np.ma.masked_array(obs.data, mask=region)
good_lvl = np.zeros(bin_num,dtype=bool) good_lvl = np.zeros(bin_num, dtype=bool)
good_idx = 0 good_idx = 0
if obs.exposure > 1000: if obs.exposure > 1000:
wav_obs = obs.wavdecomp('gauss',(5,3),occ_coeff=True) wav_obs = obs.wavdecomp('gauss', thresh, occ_coeff=True)
occ_coeff = obs.get_coeff() occ_coeff = obs.get_coeff()
for idx, lvl in enumerate(binary_array(bin_num)): for idx, lvl in enumerate(binary_array(bin_num)):
try: try:
region = wav_obs[2:-1][lvl].sum(0)>0 region = wav_obs[2:-1][lvl].sum(0) > 0
except ValueError: except ValueError:
region = np.zeros(obs.data.shape,dtype=bool) region = np.zeros(obs.data.shape, dtype=bool)
masked_obs = np.ma.masked_array(obs.data, mask = region)*occ_coeff masked_obs = np.ma.masked_array(obs.data, mask=region)*occ_coeff
rem_region = np.logical_and(region, np.logical_not(obs.data.mask)) rem_region = np.logical_and(region, np.logical_not(obs.data.mask))
rem_signal[idx] = 1-obs.data[region].sum()/obs.data.sum() rem_signal[idx] = 1-obs.data[region].sum()/obs.data.sum()
rem_area[idx] = 1 - rem_region.sum()/np.logical_not(obs.data.mask).sum() rem_area[idx] = 1 - rem_region.sum()/np.logical_not(obs.data.mask).sum()
poiss_comp[idx] = np.mean((masked_obs-masked_obs.mean())**2/masked_obs.mean()) poiss_comp[idx] = np.mean((masked_obs-masked_obs.mean())**2/masked_obs.mean())
rms[idx] = np.sqrt(((masked_obs-masked_obs.mean())**2).mean()) rms[idx] = np.sqrt(((masked_obs-masked_obs.mean())**2).mean())
parameter = lambda idx: ((poiss_comp[idx])**2+((1-rem_area[idx])*0.5)**2) parameter = lambda idx: ((poiss_comp[idx])**2+((1-rem_area[idx])*0.5)**2)
if (parameter(idx)<parameter(good_idx)): if (parameter(idx) < parameter(good_idx)):
good_idx = idx good_idx = idx
good_lvl = lvl good_lvl = lvl
try: try:
region = wav_obs[2:-1][good_lvl].sum(0)>0 region = wav_obs[2:-1][good_lvl].sum(0) > 0
except ValueError: except ValueError:
region = np.zeros(obs.data.shape,dtype=bool) region = np.zeros(obs.data.shape, dtype=bool)
masked_obs = np.ma.masked_array(obs.data, mask = region) region_raw = obs.region_to_raw(region.astype(int))
masked_obs = np.ma.masked_array(obs.data, mask=region)
rem_region = np.logical_and(region, np.logical_not(obs.data.mask)) rem_region = np.logical_and(region, np.logical_not(obs.data.mask))
to_table = [obs.obs_id, to_table = [obs.obs_id,
obs.det, obs.det,
@ -251,8 +315,8 @@ def process(obs_path):
lat, lat,
obs.time_start, obs.time_start,
obs.exposure, obs.exposure,
masked_obs.mean()/obs.exposure, #count rate masked_obs.mean()/obs.exposure, # count rate
1 - rem_region.sum()/np.logical_not(obs.data.mask).sum(), #rem_area 1 - rem_region.sum()/np.logical_not(obs.data.mask).sum(), # rem_area
poiss_comp[good_idx], poiss_comp[good_idx],
poiss_comp[0], poiss_comp[0],
rms[good_idx] rms[good_idx]
@ -267,89 +331,104 @@ def process(obs_path):
obs.time_start, obs.time_start,
obs.exposure, obs.exposure,
-1, -1,
-1, #rem_signal -1, # rem_signal
-1, #rem_area -1, # rem_area
-1, -1,
-1, -1,
-1 -1
] ]
return to_table, region.astype(int) return to_table, region.astype(int), region_raw
except TypeError: except TypeError:
return obs_path, np.zeros((360,360)) return obs_path, -1, -1
def process_folder(input_folder = None, continue_old = None, fits_folder = None):
#DIALOGUE
if not(input_folder): def process_folder(input_folder=None, start_new_file=None, fits_folder=None, thresh=None):
# DIALOGUE
if not (input_folder):
print('Enter path to the input folder') print('Enter path to the input folder')
input_folder = input() input_folder = input()
if not(continue_old): if not (start_new_file):
print('Create new file for this processing? y/n') print('Create new file for this processing? y/n')
continue_old = input() start_new_file = input()
if continue_old == 'y': if start_new_file == 'y':
start_new = True start_new = True
elif continue_old == 'n': elif start_new_file == 'n':
start_new = False start_new = False
else: else:
print('Cannot interprete input, closing script') print('Cannot interprete input, closing script')
raise SystemExit(0) raise SystemExit(0)
if not(fits_folder): if not (fits_folder):
print(f'Enter path to the output folder') print(f'Enter path to the output folder')
fits_folder = input() fits_folder = input()
region_folder = f'{fits_folder}\\Region' region_folder = f'{fits_folder}\\Region'
#CREATE ALL NECESSARY FILES AND VARIBALES region_raw_folder = f'{fits_folder}\\Region_raw'
obs_list = get_link_list(input_folder,sort_list = True) if not thresh:
print('Enter threshold values for wavelet decomposition:')
print('General threshold:')
_thresh_max = float(input())
print('Additional threshold:')
_thresh_add = float(input())
thresh = (_thresh_max, _thresh_add)
# CREATE ALL NECESSARY FILES AND VARIBALES
obs_list = get_link_list(input_folder, sort_list=True)
start = perf_counter() start = perf_counter()
group_size = 50 group_size = 50
makedirs(region_folder,exist_ok = True) makedirs(region_folder, exist_ok=True)
#FILTERING BY THE FILE SIZE makedirs(region_raw_folder, exist_ok=True)
# FILTERING BY THE FILE SIZE
print(f'Finished scanning folders. Found {len(obs_list)} observations.') print(f'Finished scanning folders. Found {len(obs_list)} observations.')
table = { table = {
'obs_id':[], 'detector':[], 'ra':[], 'dec':[], 'lon':[], 'lat':[], 't_start':[], 'exposure':[], 'obs_id': [], 'detector': [], 'ra': [], 'dec': [],
'count_rate':[], 'remaining_area':[], 'poisson_chi2':[], 'poisson_chi2_full':[], 'rms':[] 'lon': [], 'lat': [], 't_start': [], 'exposure': [],
'count_rate': [], 'remaining_area': [], 'poisson_chi2': [],
'poisson_chi2_full': [], 'rms': []
} }
if start_new: if start_new:
out_table = DataFrame(table) out_table = DataFrame(table)
out_table.to_csv(f'{fits_folder}\\test.csv') out_table.to_csv(f'{fits_folder}\\test.csv')
out_table.to_csv(f'{fits_folder}\\test_skipped.csv') out_table.to_csv(f'{fits_folder}\\test_skipped.csv')
#FILTERING OUT PROCESSED OBSERVATIONS # FILTERING OUT PROCESSED OBSERVATIONS
already_processed_list = read_csv(f'{fits_folder}\\test.csv',index_col=0,dtype={'obs_id':str}) already_processed_list = read_csv(f'{fits_folder}\\test.csv', index_col=0, dtype={'obs_id':str})
already_skipped_list = read_csv(f'{fits_folder}\\test_skipped.csv',index_col=0,dtype={'obs_id':str}) already_skipped_list = read_csv(f'{fits_folder}\\test_skipped.csv', index_col=0, dtype={'obs_id':str})
already_processed = (already_processed_list['obs_id'].astype(str)+already_processed_list['detector']).values already_processed = (already_processed_list['obs_id'].astype(str)+already_processed_list['detector']).values
already_skipped = (already_skipped_list['obs_id'].astype(str)+already_skipped_list['detector']).values already_skipped = (already_skipped_list['obs_id'].astype(str)+already_skipped_list['detector']).values
obs_list_names = [curr[curr.index('nu')+2:curr.index('_cl.evt')-2] for curr in obs_list] obs_list_names = [curr[curr.index('nu')+2:curr.index('_cl.evt')-2] for curr in obs_list]
not_processed = np.array([(curr not in already_processed) for curr in obs_list_names]) not_processed = np.array([(curr not in already_processed) for curr in obs_list_names])
not_skipped = np.array([(curr not in already_skipped) for curr in obs_list_names]) not_skipped = np.array([(curr not in already_skipped) for curr in obs_list_names])
obs_list = obs_list[np.logical_and(not_processed,not_skipped)] obs_list = obs_list[np.logical_and(not_processed, not_skipped)]
print(f'Removed already processed observations. {len(obs_list)} observations remain.') print(f'Removed already processed observations. {len(obs_list)} observations remain.')
#START PROCESSING # START PROCESSING
print('Started processing...') print('Started processing...')
num = 0 num = 0
for group_idx in range(len(obs_list)//group_size+1): for group_idx in range(len(obs_list)//group_size+1):
print(f'Started group {group_idx}') print(f'Started group {group_idx}')
group_list = obs_list[group_size*group_idx:min(group_size*(group_idx+1),len(obs_list))] group_list = obs_list[group_size*group_idx:min(group_size*(group_idx+1), len(obs_list))]
max_size = np.array([stat(file).st_size/2**20 for file in group_list]).max() max_size = np.array([stat(file).st_size/2**20 for file in group_list]).max()
process_num = cpu_count() if max_size<50 else (cpu_count()//2 if max_size<200 else (cpu_count()//4 if max_size<1000 else 1)) process_num = cpu_count() if max_size < 50 else (cpu_count()//2 if max_size < 200 else (cpu_count()//4 if max_size < 1000 else 1))
print(f"Max file size in group is {max_size:.2f}Mb, create {process_num} processes") print(f"Max file size in group is {max_size:.2f}Mb, create {process_num} processes")
with get_context('spawn').Pool(processes=process_num) as pool: with get_context('spawn').Pool(processes=process_num) as pool:
for result,region in pool.imap(process,group_list): packed_args = map(lambda _: (_, thresh), group_list)
for result, region, region_raw in pool.imap(process, packed_args):
if type(result) is np.str_: if type(result) is np.str_:
obs_id = result[result.index('nu'):result.index('_cl.evt')] obs_id = result[result.index('nu'):result.index('_cl.evt')]
print(f'{num:>3} is skipped. File {obs_id}') print(f'{num:>3} is skipped. File {obs_id}')
num +=1 num += 1
continue continue
for key,value in zip(table.keys(),result): for key, value in zip(table.keys(), result):
table[key] = [value] table[key] = [value]
if table['exposure'][0] < 1000: if table['exposure'][0] < 1000:
print(f'{num:>3} {str(result[0])+result[1]} is skipped. Exposure < 1000') print(f'{num:>3} {str(result[0])+result[1]} is skipped. Exposure < 1000')
DataFrame(table).to_csv(f'{fits_folder}\\test_skipped.csv',mode='a',header=False) DataFrame(table).to_csv(f'{fits_folder}\\test_skipped.csv', mode='a', header=False)
num +=1 num +=1
continue continue
DataFrame(table).to_csv(f'{fits_folder}\\test.csv',mode='a',header=False) DataFrame(table).to_csv(f'{fits_folder}\\test.csv', mode='a', header=False)
fits.writeto(f'{region_folder}\\{str(result[0])+result[1]}_region.fits', region, overwrite= True) fits.writeto(f'{region_folder}\\{str(result[0])+result[1]}_region.fits', region, overwrite=True)
fits.writeto(f'{region_raw_folder}\\{str(result[0])+result[1]}_reg_raw.fits', region_raw, overwrite=True)
print(f'{num:>3} {str(result[0])+result[1]} is written.') print(f'{num:>3} {str(result[0])+result[1]} is written.')
num +=1 num +=1
print('Converting generated csv to fits file...') print('Converting generated csv to fits file...')
print(f'Current time in: {(perf_counter()-start):.2f}') print(f'Current time in: {(perf_counter()-start):.2f}')
print(f'Processed {num/len(obs_list)*100:.2f} percent') print(f'Processed {num/len(obs_list)*100:.2f} percent')
csv_file = read_csv(f'{fits_folder}\\test.csv',index_col=0,dtype={'obs_id':str}) csv_file = read_csv(f'{fits_folder}\\test.csv', index_col=0, dtype={'obs_id': str})
Table.from_pandas(csv_file).write(f'{fits_folder}\\test.fits',overwrite=True) Table.from_pandas(csv_file).write(f'{fits_folder}\\test.fits', overwrite=True)
print(f'Finished writing: {perf_counter()-start}') print(f'Finished writing: {perf_counter()-start}')