From 9e2c094b0ad096489d426844eb0883be966a30f3 Mon Sep 17 00:00:00 2001 From: Andrey Mukhin Date: Wed, 31 Aug 2022 16:23:09 +0300 Subject: [PATCH 1/9] commit --- get_region_archive_table.py | 14 +---- get_region_pack.py | 101 ++++++++++++------------------------ 2 files changed, 34 insertions(+), 81 deletions(-) diff --git a/get_region_archive_table.py b/get_region_archive_table.py index 7a24949..5b9a492 100644 --- a/get_region_archive_table.py +++ b/get_region_archive_table.py @@ -12,16 +12,6 @@ import time import os warnings.filterwarnings('ignore') # %% -def ellipse(array): - grid = np.indices(array.shape)[::-1] - center = [((curr_axis*array).sum()/array.sum()) for curr_axis in grid] - y,x = np.where(array) - return center, np.abs(x-center[0]).max()*2, np.abs(y-center[1]).max()*2 -def mask_ellipse(array): - (center_x, center_y), len_x, len_y = ellipse(array) - x,y = np.indices(array.shape)[::-1] - radius = ((x-center_x)/(0.5*len_x))**2+((y-center_y)/(0.5*len_y))**2 - return radius <= 1 def poisson_divider(array): sub_arrays = np.zeros((4,180,180)) for i in range(2): @@ -109,8 +99,8 @@ def process(argument): #%% if __name__ == '__main__': start = time.perf_counter() - processing = True - start_new = True + processing = False + start_new = False group_size = 50 fits_folder = 'D:\Programms\Jupyter\Science\Source_mask\\\Archive\Processing_v8' region_folder = f'{fits_folder}\\Region' diff --git a/get_region_pack.py b/get_region_pack.py index 41243c0..a00a7f7 100644 --- a/get_region_pack.py +++ b/get_region_pack.py @@ -60,25 +60,12 @@ def atrous(level = 0, resize = False, max_size = 1001): 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 def gauss(level=0, resize=False, max_size = 1000): - size = min(5*2**level+1 if not resize else 5*2**(level+1)+1, max_size) + size = min(5*2**(level+1)+1, max_size) sigma = 2**(level) - if sigma < 1: - out = np.zeros((size+1,size+1)) - out[int((size-1)/2)+1][int((size-1)/2)+1] = 1 - return out A = 1/(2*np.pi*sigma**2)**0.5 x = A*np.exp( (-(np.arange(size)-(size-1)//2)**2)/(2*sigma**2)) out = np.multiply.outer(x,x) return out -def sigma(mode='atrous',level=0): - if mode=='atrous': - sigma = [0.8908, 0.20066, 0.08551, 0.04122, 0.02042, 0.0114] - elif mode=='gauss': - sigma = [0.912579, 0.125101, 0.104892, 4.97810e-02, 2.46556e-02, 1.14364e-02] - if level < 6: - return sigma[level] - else: - return sigma[5]/2**(level-5) def adjecent(array): grid = np.array([ [1,1,1], @@ -166,31 +153,10 @@ class Observation: self.model = (fits.getdata(f'D:\Programms\Jupyter\Science\Source_mask\Model/det1_fpm{self.det}.cxb.fits',0) *(180/np.pi/10150)**2)[shift:360+shift,shift:360+shift] self.exposure = self.header['EXPOSURE'] - # def get_coeff(self,folder='D:\\Programms\\Jupyter\\Science\\Source_mask\\Archive\\OCC_observations'): - # try: - # # _occ_list = get_link_list('D:\\Programms\\Jupyter\\Science\\Source_mask\\Archive\\OCC_observations') - # _occ_list = get_link_list(folder) - # occ_list = dict() - # for name in _occ_list: - # occ_list[name[name.index('nu'):name.index('02_cl.evt')]] = name - # occ = Observation(occ_list[self.name[:-2]],E_borders=[10,20]) - # output = np.zeros((360,360)) - # for i in range(2): - # for j in range(2): - # _temp = occ.data[180*i:180*(i+1),180*j:180*(j+1)].mean() - # output[180*i:180*(i+1),180*j:180*(j+1)] = _temp if _temp != 0 else 1 - # except KeyError: - # output = np.ones((360,360)) - # out = output.min()/output - # if np.any(out<0.8): - # return np.ones((360,360)) - # else: - # return out def get_coeff(self): - # coeff = np.array([0.972,0.895,1.16,1.02]) if self.det=='A' else np.array([0.991,1.047,1.012,0.956]) 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]) - test = (coeff).reshape(2,2).repeat(180,0).repeat(180,1) - return test + resized_coeff = (coeff).reshape(2,2).repeat(180,0).repeat(180,1) + return resized_coeff def get_data(self, file, E_borders=[3,20]): PI_min, PI_max = (np.array(E_borders)-1.6)/0.04 data = file[1].data.copy() @@ -215,15 +181,14 @@ class Observation: temp = np.logical_or(temp, np.equal(pixpos_file['rawx'],x)*np.equal(pixpos_file['rawy'],y)) temp = pixpos_file[temp] 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 - def wavdecomp(self, mode = 'atrous', thresh=False, iteration = 1,occ_coeff = False): + def wavdecomp(self, mode = 'atrous', thresh=False,occ_coeff = False): #THRESHOLD if type(thresh) is int: thresh_max, thresh_add = thresh, thresh/2 elif type(thresh) is tuple: thresh_max, thresh_add = thresh - #INIT BASIC VARIABLES + #INIT NEEDED VARIABLES wavelet = globals()[mode] - # max_level = int(np.ceil(np.log2(self.data.shape[0])))+1 max_level = 8 conv_out = np.zeros((max_level+1,self.data.shape[0],self.data.shape[1])) size = self.data.shape[0] @@ -234,36 +199,34 @@ class Observation: data_bkg = data.copy() #ITERATIVELY CONDUCT WAVLET DECOMPOSITION for i in range(max_level): - for num_iter in range(iteration): - conv = fftconvolve(data,wavelet(i),mode='same') - temp_out = data-conv - #ERRORMAP CALCULATION - if thresh_max != 0: - sig = sigma(mode, i) - bkg = fftconvolve(data_bkg, wavelet(i),mode='same') - bkg[bkg<0] = 0 - # err = (1+np.sqrt(bkg/sig**2 + 0.75))*sig**3 - 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] - 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] + conv = fftconvolve(data,wavelet(i),mode='same') + temp_out = data-conv + #ERRORMAP CALCULATION + if thresh_max != 0: + sig = sigma(mode, i) + bkg = fftconvolve(data_bkg, wavelet(i),mode='same') + bkg[bkg<0] = 0 + # err = (1+np.sqrt(bkg/sig**2 + 0.75))*sig**3 + 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] + 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] + adj = adjecent(significant) + add_condition = np.logical_and(add_significant[adj[0],adj[1]],np.logical_not(significant[adj[0],adj[1]])) + while (add_condition).any(): + to_add = adj[0][add_condition], adj[1][add_condition] + significant[to_add[0],to_add[1]] = True adj = adjecent(significant) add_condition = np.logical_and(add_significant[adj[0],adj[1]],np.logical_not(significant[adj[0],adj[1]])) - while (add_condition).any(): - to_add = adj[0][add_condition], adj[1][add_condition] - significant[to_add[0],to_add[1]] = True - 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(np.abs(temp_out)[adj[0],adj[1]] >= thresh_add*err[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 - #WRITING THE WAVELET DECOMP LAYER - if temp_out[size:2*size,size:2*size].sum() == 0: break - 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 - # conv_out[i] = significant - data = conv + # 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]]) + temp_out[size:2*size,size:2*size][np.logical_not(significant)] = 0 + #WRITING THE WAVELET DECOMP LAYER + if temp_out[size:2*size,size:2*size].sum() == 0: break + 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 + data = conv conv_out[max_level] = conv[size:2*size,size:2*size] return conv_out # %% From 486d655235c4713f1d775bc423f66fcee6ec519a Mon Sep 17 00:00:00 2001 From: Andrey Mukhin Date: Wed, 31 Aug 2022 16:23:44 +0300 Subject: [PATCH 2/9] commit --- get_region_pack.py | 18 ++++++------------ 1 file changed, 6 insertions(+), 12 deletions(-) diff --git a/get_region_pack.py b/get_region_pack.py index a00a7f7..1df9b70 100644 --- a/get_region_pack.py +++ b/get_region_pack.py @@ -1,16 +1,11 @@ # %% -import import_ipynb import numpy as np -import pandas as pd import itertools -from os import listdir, mkdir, stat +from os import stat from scipy.signal import fftconvolve, convolve2d -import matplotlib.pyplot as plt -from matplotlib.colors import SymLogNorm as lognorm - from astropy.io import fits from astropy.wcs import WCS @@ -203,16 +198,16 @@ class Observation: temp_out = data-conv #ERRORMAP CALCULATION if thresh_max != 0: - sig = sigma(mode, i) + sig = ((wavelet(i)**2).sum())**0.5 bkg = fftconvolve(data_bkg, wavelet(i),mode='same') bkg[bkg<0] = 0 # err = (1+np.sqrt(bkg/sig**2 + 0.75))*sig**3 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 = (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] 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 = (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] adj = adjecent(significant) add_condition = np.logical_and(add_significant[adj[0],adj[1]],np.logical_not(significant[adj[0],adj[1]])) while (add_condition).any(): @@ -223,7 +218,6 @@ class Observation: # 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]]) temp_out[size:2*size,size:2*size][np.logical_not(significant)] = 0 #WRITING THE WAVELET DECOMP LAYER - if temp_out[size:2*size,size:2*size].sum() == 0: break 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 data = conv From 9fb09e3afa27c9bfd3b4328357451cb0fac3d6de Mon Sep 17 00:00:00 2001 From: Andrey Mukhin Date: Wed, 31 Aug 2022 16:26:20 +0300 Subject: [PATCH 3/9] Create README.md --- README.md | 1 + 1 file changed, 1 insertion(+) create mode 100644 README.md diff --git a/README.md b/README.md new file mode 100644 index 0000000..4cfd582 --- /dev/null +++ b/README.md @@ -0,0 +1 @@ +# wavsource_nustar From 588b89cc2cb6bcd503ac82ddb5ec82db3d407539 Mon Sep 17 00:00:00 2001 From: Andrey Mukhin Date: Wed, 31 Aug 2022 16:26:41 +0300 Subject: [PATCH 4/9] Create test.py --- test/test.py | 1 + 1 file changed, 1 insertion(+) create mode 100644 test/test.py diff --git a/test/test.py b/test/test.py new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/test/test.py @@ -0,0 +1 @@ + From 3893d4151e351e2a1edd671c8ecb5b1b10b71a07 Mon Sep 17 00:00:00 2001 From: Andrey Mukhin Date: Wed, 31 Aug 2022 16:34:21 +0300 Subject: [PATCH 5/9] Delete test directory --- test/test.py | 1 - 1 file changed, 1 deletion(-) delete mode 100644 test/test.py diff --git a/test/test.py b/test/test.py deleted file mode 100644 index 8b13789..0000000 --- a/test/test.py +++ /dev/null @@ -1 +0,0 @@ - From aef4d43e34468b1cae27ab4dce68d3f1dbb469ed Mon Sep 17 00:00:00 2001 From: Andrey Mukhin Date: Wed, 31 Aug 2022 16:35:22 +0300 Subject: [PATCH 6/9] new file: test.py modified: test.py --- test.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 test.py diff --git a/test.py b/test.py new file mode 100644 index 0000000..e69de29 From 5ebcc8bafa36919c476009de8b571504a1c22c1a Mon Sep 17 00:00:00 2001 From: Andrey Mukhin Date: Wed, 31 Aug 2022 16:36:38 +0300 Subject: [PATCH 7/9] commit --- test.py | 1 + 1 file changed, 1 insertion(+) diff --git a/test.py b/test.py index e69de29..7ddfc8f 100644 --- a/test.py +++ b/test.py @@ -0,0 +1 @@ +print('hello world!') \ No newline at end of file From c6f6c3402614f217f94a0fd9e1278fa56e05de32 Mon Sep 17 00:00:00 2001 From: Andrey Mukhin Date: Wed, 31 Aug 2022 16:38:14 +0300 Subject: [PATCH 8/9] Delete test.py --- test.py | 1 - 1 file changed, 1 deletion(-) delete mode 100644 test.py diff --git a/test.py b/test.py deleted file mode 100644 index 7ddfc8f..0000000 --- a/test.py +++ /dev/null @@ -1 +0,0 @@ -print('hello world!') \ No newline at end of file From fc24392f8f8dc21f7424463152f171ad8a2efa0b Mon Sep 17 00:00:00 2001 From: Andrey Mukhin Date: Wed, 31 Aug 2022 16:46:16 +0300 Subject: [PATCH 9/9] commit --- test.py | 1 - 1 file changed, 1 deletion(-) delete mode 100644 test.py diff --git a/test.py b/test.py deleted file mode 100644 index 7ddfc8f..0000000 --- a/test.py +++ /dev/null @@ -1 +0,0 @@ -print('hello world!') \ No newline at end of file