diff --git a/box_helpers/move_boxes.py b/box_helpers/move_boxes.py index 065e2790..e942e7b6 100644 --- a/box_helpers/move_boxes.py +++ b/box_helpers/move_boxes.py @@ -6,6 +6,8 @@ from glob import glob import os +from argparse import ArgumentParser + def move_boxes(file, folder): try: @@ -30,8 +32,9 @@ def move_boxes(file, folder): print("Failing to open ds9 to verify box selection, check if installed and try to run on the commandline" "\nds9 {FILE} -regions load all '{DATALOC}/boxes/*.reg'".format(FILE=file, DATALOC=folder)) -if __name__ == '__main__': - from argparse import ArgumentParser +def parse_args(): + """Command line parser""" + parser = ArgumentParser() parser.add_argument('-l', '--location', type=str, help='data location folder name', default='.') parser.add_argument('-f', '--file', type=str, help='fitsfile name', @@ -42,4 +45,14 @@ def move_boxes(file, folder): if folder[-1] == '/': folder = folder[0:-1] + return args, folder + +def main(): + """Main function""" + + args, folder = parse_args() + move_boxes(args.file, folder) + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/h5_helpers/add_h5_dirs.py b/h5_helpers/add_h5_dirs.py index 7bbcac72..d059f856 100644 --- a/h5_helpers/add_h5_dirs.py +++ b/h5_helpers/add_h5_dirs.py @@ -3,17 +3,30 @@ import os import numpy as np -if __name__ == '__main__': + +def parse_args(): + """ + Command line argument parser + + :return: parsed arguments + """ parser = ArgumentParser() - parser.add_argument('--h5_in', type=str, help='input h5 (from which to extract the frequencies and antennas)', required=True) + parser.add_argument('--h5_in', type=str, help='input h5 (from which to extract the frequencies and antennas)', + required=True) parser.add_argument('--h5_dirs', type=str, help='h5 with directions to copy from', required=True) parser.add_argument('--h5_out', type=str, help='output h5') - args = parser.parse_args() + return parser.parse_args() + + +def main(): + """Main function""" + + args = parse_args() os.system(f'cp {args.h5_in} {args.h5_out}') - print("Make "+args.h5_out) + print("Make " + args.h5_out) H = tables.open_file(args.h5_out, "r+") T = tables.open_file(args.h5_dirs) @@ -41,4 +54,8 @@ H.root.sol000._f_get_child(table)._f_get_child(axes).attrs['AXES'] = AXES H.root.sol000._f_get_child(table).dir._f_remove() - H.create_array(H.root.sol000._f_get_child(table), 'dir', T.root.sol000._f_get_child(table).dir[:]) \ No newline at end of file + H.create_array(H.root.sol000._f_get_child(table), 'dir', T.root.sol000._f_get_child(table).dir[:]) + + +if __name__ == '__main__': + main() diff --git a/h5_helpers/find_closest_h5.py b/h5_helpers/find_closest_h5.py index 20d79f28..e9e8cc10 100644 --- a/h5_helpers/find_closest_h5.py +++ b/h5_helpers/find_closest_h5.py @@ -5,10 +5,12 @@ import sys from argparse import ArgumentParser + class FindClosestDir: """ Make template h5 """ + def __init__(self, h5_in, template_name): os.system(' '.join(['cp', h5_in, template_name])) print(f'Created {template_name}') @@ -20,7 +22,6 @@ def __init__(self, h5_in, template_name): self.dirs = T.root.sol000.source[:]['dir'] T.close() - def make_template(self): """ Make template h5 with 1 direction @@ -32,7 +33,7 @@ def make_template(self): for solset in self.h5.root._v_groups.keys(): ss = self.h5.root._f_get_child(solset) ss._f_get_child('source')._f_remove() - values = np.array([(b'Dir00', [0. , 0. ])], dtype=[('name', 'S128'), ('dir', ' output_h5s\n' 'You can optionally use h5_merger.py to merge all the directions into 1 big h5 file\n' 'Example command --> python h5_merger.py -in output_h5s/source_*.h5 -out merged.h5 --propagate_flags') + + +if __name__ == '__main__': + main() diff --git a/h5_helpers/h5_filter.py b/h5_helpers/h5_filter.py index 8e4ca8da..324fbd36 100644 --- a/h5_helpers/h5_filter.py +++ b/h5_helpers/h5_filter.py @@ -16,6 +16,7 @@ from glob import glob import re + def str2bool(v): v = str(v) if v.lower() in ('yes', 'true', 't', 'y', '1'): @@ -25,24 +26,29 @@ def str2bool(v): else: raise ArgumentTypeError('Boolean value expected.') + def degree_to_radian(inp): """degree to radian""" return float(inp) / 360 * pi * 2 + def radian_to_degree(inp): """radion to degree""" return float(inp) * 360 / (pi * 2) + def angular_distance(p1, p2): """angular distance for points in ra and dec in degrees""" - if p1[0]>2*pi: + if p1[0] > 2 * pi: p1 = [degree_to_radian(p) for p in p1] p2 = [degree_to_radian(p) for p in p2] - return radian_to_degree(acos(sin(p1[1])*sin(p2[1])+cos(p1[1])*cos(p2[1])*cos(p1[0]-p2[0]))) + return radian_to_degree(acos(sin(p1[1]) * sin(p2[1]) + cos(p1[1]) * cos(p2[1]) * cos(p1[0] - p2[0]))) + def remove_numbers(inp): return "".join(re.findall("[a-zA-z]+", inp)) + def create_new_soltab(h5_in_name, h5_out_name, directions, sources): """ Create a new dataset in the h5 table @@ -68,7 +74,9 @@ def create_new_soltab(h5_in_name, h5_out_name, directions, sources): solsetout.obj.source.append(new_sources) for st in h5_in.getSolset(ss).getSoltabNames(): - print('Filter {solset}/{soltab} from {h5_in} into {h5_out}'.format(solset=ss, soltab=st, h5_in=h5_in_name.split('/')[-1], h5_out=h5_out_name.split('/')[-1])) + print('Filter {solset}/{soltab} from {h5_in} into {h5_out}'.format(solset=ss, soltab=st, + h5_in=h5_in_name.split('/')[-1], + h5_out=h5_out_name.split('/')[-1])) solutiontable = h5_in.getSolset(ss).getSoltab(st) axes = solutiontable.getValues()[1] @@ -77,65 +85,86 @@ def create_new_soltab(h5_in_name, h5_out_name, directions, sources): axes['dir'] = [ns[0] for ns in new_sources] dir_index = solutiontable.getAxesNames().index('dir') shape = list(values_in.shape) - shape[dir_index]=len(directions) + shape[dir_index] = len(directions) values_new = zeros(shape) for idx_new, idx_old in enumerate(indexes): if dir_index == 0: - values_new[idx_new,...] += values_in[idx_old, ...] + values_new[idx_new, ...] += values_in[idx_old, ...] elif dir_index == 1: - values_new[:, idx_new,...] += values_in[:, idx_old, ...] + values_new[:, idx_new, ...] += values_in[:, idx_old, ...] elif dir_index == 2: - values_new[:, :, idx_new,...] += values_in[:, :, idx_old, ...] + values_new[:, :, idx_new, ...] += values_in[:, :, idx_old, ...] elif dir_index == 3: - values_new[:, :, :, idx_new,...] += values_in[:, :, :, idx_old, ...] + values_new[:, :, :, idx_new, ...] += values_in[:, :, :, idx_old, ...] elif dir_index == 4: - values_new[:, :, :, :, idx_new,...] += values_in[:, :, :, :, idx_old,...] + values_new[:, :, :, :, idx_new, ...] += values_in[:, :, :, :, idx_old, ...] print('New number of sources {num}'.format(num=len(sources))) print('Filtered output shape {shape}'.format(shape=values_new.shape)) weights = ones(values_new.shape) - solsetout.makeSoltab(remove_numbers(st), axesNames=list(axes.keys()), axesVals=list(axes.values()), vals=values_new, - weights=weights) + solsetout.makeSoltab(remove_numbers(st), axesNames=list(axes.keys()), axesVals=list(axes.values()), + vals=values_new, + weights=weights) h5_in.close() h5_out.close() -parser = ArgumentParser() -parser.add_argument('-f', '--fits', type=str, help='fitsfile name') -parser.add_argument('-ac', '--angular_cutoff', type=float, default=None, help='angular distances higher than this value from the center will be excluded from the box selection') -parser.add_argument('-in', '--inside', type=str2bool, default=False, help='keep directions inside the angular cutoff') -parser.add_argument('-h5out', '--h5_file_out', type=str, help='h5 output name') -parser.add_argument('-h5in', '--h5_file_in', type=str, help='h5 input name (to filter)') -args = parser.parse_args() - -# clean up files -if args.h5_file_out.split('/')[-1] in [f.split('/')[-1] for f in glob(args.h5_file_out)]: - os.system('rm -rf {file}'.format(file=args.h5_file_out)) - -# get header from fits file -hdu = fits.open(args.fits)[0] -header = WCS(hdu.header, naxis=2).to_header() -# get center of field -center = (degree_to_radian(header['CRVAL1']), degree_to_radian(header['CRVAL2'])) - -# return list of directions that have to be included -H = tables.open_file(args.h5_file_in) -sources = [] -directions = [] -for dir in H.root.sol000.source[:]: - print(angular_distance(center, dir[1])) - if args.inside and angular_distance(center, dir[1])=args.angular_cutoff: - print('Keep {dir}'.format(dir=dir)) - directions.append(dir[0]) - sources.append(dir) - else: - print('Remove {dir}'.format(dir=dir)) -H.close() -create_new_soltab(args.h5_file_in, args.h5_file_out, directions, sources) \ No newline at end of file +def parse_args(): + """ + Command line argument parser + + :return: parsed arguments + """ + + parser = ArgumentParser() + parser.add_argument('-f', '--fits', type=str, help='fitsfile name') + parser.add_argument('-ac', '--angular_cutoff', type=float, default=None, + help='angular distances higher than this value from the center will be excluded from the box selection') + parser.add_argument('-in', '--inside', type=str2bool, default=False, + help='keep directions inside the angular cutoff') + parser.add_argument('-h5out', '--h5_file_out', type=str, help='h5 output name') + parser.add_argument('-h5in', '--h5_file_in', type=str, help='h5 input name (to filter)') + return parser.parse_args() + + +def main(): + """Main function""" + + args = parse_args() + + # clean up files + if args.h5_file_out.split('/')[-1] in [f.split('/')[-1] for f in glob(args.h5_file_out)]: + os.system('rm -rf {file}'.format(file=args.h5_file_out)) + + # get header from fits file + hdu = fits.open(args.fits)[0] + header = WCS(hdu.header, naxis=2).to_header() + # get center of field + center = (degree_to_radian(header['CRVAL1']), degree_to_radian(header['CRVAL2'])) + + # return list of directions that have to be included + H = tables.open_file(args.h5_file_in) + sources = [] + directions = [] + for dir in H.root.sol000.source[:]: + print(angular_distance(center, dir[1])) + if args.inside and angular_distance(center, dir[1]) < args.angular_cutoff: + print('Keep {dir}'.format(dir=dir)) + directions.append(dir[0]) + sources.append(dir) + elif not args.inside and angular_distance(center, dir[1]) >= args.angular_cutoff: + print('Keep {dir}'.format(dir=dir)) + directions.append(dir[0]) + sources.append(dir) + else: + print('Remove {dir}'.format(dir=dir)) + H.close() + + create_new_soltab(args.h5_file_in, args.h5_file_out, directions, sources) + + +if __name__ == '__main__': + main() diff --git a/h5_helpers/h5_flagger.py b/h5_helpers/h5_flagger.py index a4d6e0b6..1c7e42fa 100644 --- a/h5_helpers/h5_flagger.py +++ b/h5_helpers/h5_flagger.py @@ -3,17 +3,31 @@ from argparse import ArgumentParser import re -if __name__ == '__main__': + +def parse_args(): + """ + Command line argument parser + + :return: parsed arguments + """ parser = ArgumentParser() parser.add_argument('--h5', nargs='+', help='Input h5parm', required=True) parser.add_argument('--freqrange', help='Flag frequency range in MHz (example: 140-158') parser.add_argument('--ant', help='Antenna to flag (example: RS409HBA)') parser.add_argument('--ampflag', help='Flag based on amplitude values', action='store_true') - args = parser.parse_args() + return parser.parse_args() + + +def main(): + """ + Main script + """ + + args = parse_args() if args.freqrange is not None: - min_freq, max_freq = [float(f)*1e6 for f in args.freqrange.split('-')] + min_freq, max_freq = [float(f) * 1e6 for f in args.freqrange.split('-')] ant = args.ant print(f'Flagging from {min_freq}Hz to {max_freq}Hz for antenna {ant}') @@ -21,7 +35,7 @@ ampflag = args.ampflag for h5 in args.h5: - print('-----\n'+h5+'\n-----') + print('-----\n' + h5 + '\n-----') H = tables.open_file(h5, 'r+') for solset in H.root._v_groups.keys(): print(solset) @@ -42,24 +56,25 @@ del shape[AXES.index('ant')] newweights = np.zeros(shape) - if AXES.index('ant')==2 and AXES.index('freq')==1: + if AXES.index('ant') == 2 and AXES.index('freq') == 1: print(weights[:, freq_indices[0]:freq_indices[-1], ant_idx, ...].shape, newweights.shape) - H.root._f_get_child(solset)._f_get_child(soltab).\ - weight[:, freq_indices[0]:freq_indices[-1], ant_idx, ...] = newweights + H.root._f_get_child(solset)._f_get_child(soltab). \ + weight[:, freq_indices[0]:freq_indices[-1], ant_idx, ...] = newweights if args.ampflag is not None and 'amplitude' in soltab: values = st.val[:] flag_vals = (values < 3) & (values > 0.33) - flag_percentage = 1-(flag_vals*weights).sum()/weights.sum() - print(f'Flagged {round(flag_percentage*100, 4)}%') + flag_percentage = 1 - (flag_vals * weights).sum() / weights.sum() + print(f'Flagged {round(flag_percentage * 100, 4)}%') del values ss._f_get_child(soltab).weight[:] *= flag_vals phasenum = re.findall(r'\d+', soltab)[0] - if 'phase'+phasenum in list(ss._v_groups.keys()): - ss._f_get_child('phase'+phasenum).weight[:] *= flag_vals - + if 'phase' + phasenum in list(ss._v_groups.keys()): + ss._f_get_child('phase' + phasenum).weight[:] *= flag_vals H.close() +if __name__ == '__main__': + main() diff --git a/h5_helpers/pol_phase_rot.py b/h5_helpers/pol_phase_rot.py index 06957fcc..654db047 100644 --- a/h5_helpers/pol_phase_rot.py +++ b/h5_helpers/pol_phase_rot.py @@ -13,12 +13,14 @@ ----------------------------- """ + class PhaseRotate: """ Make template h5 with default values (phase=0 and amplitude=1) and convert to phase rotate matrix for polarization alignment between different observations. """ - def __init__(self, name_in, name_out, freqs = None): + + def __init__(self, name_in, name_out, freqs=None): os.system(' '.join(['cp', name_in, name_out])) self.h5 = tables.open_file(name_out, 'r+') self.axes = ['time', 'freq', 'ant', 'dir', 'pol'] @@ -58,12 +60,12 @@ def update_array(self, st, new_val, arrayname): self.h5.create_array(st, arrayname, new_val.astype(valtype), atom=atomtype) else: self.h5.create_array(st, arrayname, new_val.astype(valtype)) - if arrayname=='val' or arrayname=='weight': + if arrayname == 'val' or arrayname == 'weight': st._f_get_child(arrayname).attrs['AXES'] = bytes(','.join(self.axes), 'utf-8') return self - def make_template(self, shape = None, polrot = None): + def make_template(self, shape=None, polrot=None): """ Make template h5 @@ -72,7 +74,7 @@ def make_template(self, shape = None, polrot = None): """ print("1) MAKE TEMPLATE H5PARM") - amplitudedone=False + amplitudedone = False for solset in self.h5.root._v_groups.keys(): ss = self.h5.root._f_get_child(solset) for soltab in ss._v_groups.keys(): @@ -89,7 +91,7 @@ def make_template(self, shape = None, polrot = None): new_val = np.ones(shape) new_val[..., 1] = 0 new_val[..., 2] = 0 - amplitudedone=True + amplitudedone = True else: continue self.time = np.array([st.time[:][0]]) @@ -101,11 +103,10 @@ def make_template(self, shape = None, polrot = None): if not amplitudedone: pass - #TODO: add + # TODO: add return self - def circ2lin(self): """ Convert circular polarization to linear polarization @@ -120,13 +121,11 @@ def circ2lin(self): print('\n3) CONVERTING CIRCULAR TO LINEAR POLARIZATION\n' '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%' - +circ2lin_math+ + + circ2lin_math + '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') - G = self.h5.root.sol000.amplitude000.val[:] * np.exp(self.h5.root.sol000.phase000.val[:] * 1j) - XX = (G[..., 0] + G[..., -1]) XY = 1j * (G[..., 0] - G[..., -1]) YX = 1j * (G[..., -1] - G[..., 0]) @@ -171,12 +170,12 @@ def rotate(self, intercept, rotation_measure): """ print('\n2) ADD PHASE ROTATION') - phaserot = intercept+rotation_measure*(speed_of_light/self.freqs)**2 + phaserot = intercept + rotation_measure * (speed_of_light / self.freqs) ** 2 mapping = list(zip(list(self.freqs), list(phaserot))) print('########################\nFrequency to rotation in radian (circular base):\n------------------------') for element in mapping: - print(str(int(element[0]))+'Hz --> '+str(round(element[1], 3))+'rad') + print(str(int(element[0])) + 'Hz --> ' + str(round(element[1], 3)) + 'rad') # print('Rotate with rotation angle: '+str(intercept) + ' radian') for solset in self.h5.root._v_groups.keys(): @@ -194,15 +193,24 @@ def rotate(self, intercept, rotation_measure): return self -if __name__ == '__main__': +def parse_args(): + """ + Command line argument parser + + :return: parsed arguments + """ parser = ArgumentParser() - parser.add_argument('--h5_in', type=str, help='input h5 (from which to extract the frequencies and antennas)', required=True) + parser.add_argument('--h5_in', type=str, help='input h5 (from which to extract the frequencies and antennas)', + required=True) parser.add_argument('--h5_out', type=str, help='output name (output solution file)', required=True) parser.add_argument('--intercept', type=float, help='intercept for rotation angle') parser.add_argument('--RM', type=float, help='rotation measure') - args = parser.parse_args() + return parser.parse_args() + +def main(): + args = parse_args() phaserot = PhaseRotate(args.h5_in, args.h5_out) if args.intercept is None: phaserot.make_template() @@ -211,3 +219,7 @@ def rotate(self, intercept, rotation_measure): phaserot.rotate(intercept=args.intercept, rotation_measure=args.RM) phaserot.h5.close() + + +if __name__ == '__main__': + main() diff --git a/h5_helpers/selfcal_quality.py b/h5_helpers/selfcal_quality.py new file mode 100644 index 00000000..894e7f07 --- /dev/null +++ b/h5_helpers/selfcal_quality.py @@ -0,0 +1,442 @@ +""" +This script is a supporting script to quantify the quality of the self-calibration output from facetselfcal.py from +https://github.com/rvweeren/lofar_facet_selfcal/blob/main/facetselfcal.py + +It will return plots quantifying the change of solutions and images over self-calibration cycles. +""" + +__author__ = "Jurjen de Jong (jurjendejong@strw.leidenuniv.nl)" + +import re +import tables +from glob import glob +from scipy.stats import circstd, linregress, circmean +import numpy as np +import matplotlib.pyplot as plt +from argparse import ArgumentParser +from astropy.io import fits +import csv +from skimage.filters.rank import entropy +from skimage.morphology import disk +import pandas as pd + + +class SelfcalQuality: + def __init__(self, folder: str = None, remote_only: bool = False, international_only: bool = False): + """ + Determine quality of selfcal from facetselfcal.py + + :param folder: path to where selfcal ran + :param remote_only: consider remote stations only for amp/phase stability + :param international_only: consider international stations only for amp/phase stability + """ + # selfcal folder + self.folder = folder + # merged selfcal h5parms + self.h5s = glob(f"{self.folder}/merged_selfcalcyle*.h5") + assert len(self.h5s) != 0, "No h5 files found" + + # select all sources + regex = "merged_selfcalcyle\d{3}\_" + self.sources = set([re.sub(regex, '', h.split('/')[-1]).replace('.ms.copy.phaseup.h5', '') for h in self.h5s]) + self.sourcename = list(self.sources)[0].split('_')[-1] + + # select all fits images + fitsfiles = sorted(glob(self.folder + "/*MFS-I-image.fits")) + if len(fitsfiles) == 0: + fitsfiles = sorted(glob(self.folder + "/*MFS-image.fits")) + self.fitsfiles = [f for f in fitsfiles if 'arcsectaper' not in f] + + # select all fits model images + modelfiles = sorted(glob(self.folder + "/*MFS-I-model.fits")) + if len(modelfiles) == 0: + modelfiles = sorted(glob(self.folder + "/*MFS-model.fits")) + self.modelfiles = [f for f in modelfiles if 'arcsectaper' not in f] + + # for phase/amp evolution + self.remote_only = remote_only + self.international_only = international_only + + self.textfile = open('selfcal_performance.csv', 'w') + self.writer = csv.writer(self.textfile) + self.writer.writerow(['solutions', 'dirty'] + [str(i) for i in range(len(self.fitsfiles))]) + + @staticmethod + def get_cycle_num(fitsfile: str = None): + """ + Parse cycle number + + :param fitsfile: fits file name + """ + return int(float(re.findall("\d{3}", fitsfile.split('/')[-1])[0])) + + @staticmethod + def make_utf8(inp=None): + """ + Convert input to utf8 instead of bytes + + :param inp: string input + """ + try: + inp = inp.decode('utf8') + return inp + except (UnicodeDecodeError, AttributeError): + return inp + + @staticmethod + def make_figure(vals1=None, vals2=None, label1=None, label2=None, plotname=None): + """ + Make figure (with optionally two axis) + + :param phase_scores: + :param amp_scores: + :param plotname: + :return: + """ + fig, ax1 = plt.subplots() + + color = 'tab:red' + ax1.set_xlabel('cycle') + ax1.set_ylabel(label1, color=color) + ax1.plot([i + 1 for i in range(len(vals1))], vals1, color=color) + ax1.tick_params(axis='y', labelcolor=color) + # ax1.set_ylim(0, np.pi/2) + ax1.set_xlim(1, 11) + ax1.grid(False) + ax1.grid('off') + ax1.grid(None) + + if vals2 is not None: + + ax2 = ax1.twinx() + + color = 'tab:blue' + ax2.set_ylabel(label2, color=color) + if 'Amp' in label2: + ax2.plot([i for i in range(len(vals2))][3:], vals2[3:], color=color) + else: + ax2.plot([i for i in range(len(vals2))], vals2, color=color) + ax2.tick_params(axis='y', labelcolor=color) + # ax2.set_ylim(0, 2) + ax2.set_xlim(1, 11) + ax2.grid(False) + ax2.grid('off') + ax2.grid(None) + + fig.tight_layout() + + plt.savefig(plotname, dpi=300) + + def linreg_slope(self, values=None): + """ + Fit linear regression and return slope + + :param values: Values + + :return: linear regression slope + """ + return linregress(list(range(len(values))), values).slope + + @staticmethod + def image_entropy(fitsfile: str = None): + """ + Calculate entropy of image + + :param fitsfile: + + :return: image entropy value + """ + file = fits.open(fitsfile) + image = file[0].data + file.close() + while image.ndim > 2: + image = image[0] + val = entropy(image, disk(10)).sum() / image.max() + print(val) + return val + + @staticmethod + def euclidean_distance(l1=None, l2=None): + """ + Take euclidean distance + + :return: euclidean distance + """ + return np.sqrt(np.sum(np.power(np.subtract(l1, l2), 2))) + + def get_solution_scores(self, h5_1: str = None, h5_2: str = None): + """ + Get solution scores + + :param h5_1: solution file 1 + :param h5_2: solution file 2 + + :return: phasescore --> circular std phase difference score + ampscore --> std amp difference score + """ + + # PHASE VALUES + H = tables.open_file(h5_1) + axes = self.make_utf8(H.root.sol000.phase000.val.attrs['AXES']).split(',') + vals1 = H.root.sol000.phase000.val[:] + weights1 = H.root.sol000.phase000.weight[:] + + if h5_2 is not None: + F = tables.open_file(h5_2) + vals2 = F.root.sol000.phase000.val[:] + weights2 = F.root.sol000.phase000.weight[:] + + if self.remote_only and not self.international_only: + stations = [i for i, station in enumerate(H.root.sol000.antenna[:]['name']) if + ('RS' in self.make_utf8(station))] + vals1 = np.take(vals1, stations, axis=axes.index('ant')) + weights1 = np.take(weights1, stations, axis=axes.index('ant')) + if h5_2 is not None: + vals2 = np.take(vals2, stations, axis=axes.index('ant')) + weights2 = np.take(weights2, stations, axis=axes.index('ant')) + + elif self.international_only and not self.remote_only: + stations = [i for i, station in enumerate(H.root.sol000.antenna[:]['name']) if + not ('RS' in self.make_utf8(station) + or 'CS' in self.make_utf8(station) + or 'ST' in self.make_utf8(station))] + vals1 = np.take(vals1, stations, axis=axes.index('ant')) + weights1 = np.take(weights1, stations, axis=axes.index('ant')) + if h5_2 is not None: + vals2 = np.take(vals2, stations, axis=axes.index('ant')) + weights2 = np.take(weights2, stations, axis=axes.index('ant')) + + # take circular std from difference of previous and current selfcal cycle + if h5_2 is not None: + prepphasescore = np.subtract(np.nan_to_num(vals1) * weights1, np.nan_to_num(vals2) * weights2) + else: + prepphasescore = np.nan_to_num(vals1) * weights1 + phasescore = circstd(prepphasescore[prepphasescore != 0], nan_policy='omit') + + # AMP VALUES + vals1 = H.root.sol000.amplitude000.val[:] + if h5_2 is not None: + vals2 = F.root.sol000.amplitude000.val[:] + + if self.remote_only and not self.international_only: + vals1 = np.take(vals1, stations, axis=axes.index('ant')) + if h5_2 is not None: + vals2 = np.take(vals2, stations, axis=axes.index('ant')) + + elif self.international_only and not self.remote_only: + vals1 = np.take(vals1, stations, axis=axes.index('ant')) + if h5_2 is not None: + vals2 = np.take(vals2, stations, axis=axes.index('ant')) + + # take std from ratio of previous and current selfcal cycle + if h5_2 is not None: + prepampscore = np.nan_to_num(np.divide(np.nan_to_num(vals1) * weights1, np.nan_to_num(vals2) * weights2), + posinf=0, neginf=0) + else: + prepampscore = np.nan_to_num(vals1) * weights1 + ampscore = np.std(prepampscore[prepampscore != 0]) + + H.close() + if h5_2 is not None: + F.close() + + return phasescore, ampscore + + def solution_stability(self): + """ + Get solution stability scores and make figure + + :return: bestcycle --> best cycle according to solutions + accept --> accept this selfcal + """ + + # loop over sources to get scores + for k, source in enumerate(self.sources): + print(source) + sub_h5s = sorted([h5 for h5 in self.h5s if source in h5]) + phase_scores = [] + amp_scores = [] + for m, sub_h5 in enumerate(sub_h5s): + number = self.get_cycle_num(sub_h5) + # print(sub_h5, sub_h5s[m - 1]) + if number > 0: + phasescore, ampscore = self.get_solution_scores(sub_h5, sub_h5s[m - 1]) + elif number == 0: + phasescore, ampscore = self.get_solution_scores(sub_h5, None) + phase_scores.append(phasescore) + amp_scores.append(ampscore) + if k == 0: + total_phase_scores = [phase_scores] + total_amp_scores = [amp_scores] + + total_phase_scores = np.append(total_phase_scores, [phase_scores], axis=0) + total_amp_scores = np.append(total_amp_scores, [amp_scores], axis=0) + + # plot + if self.remote_only: + plotname = f'selfcal_stability_remote_{self.sourcename}.png' + elif self.international_only: + plotname = f'selfcal_stability_international_{self.sourcename}.png' + else: + plotname = f'selfcal_stability_{self.sourcename}.png' + finalphase = np.mean(total_phase_scores, axis=0) + finalamp = np.mean(total_amp_scores, axis=0) + + self.make_figure(finalphase, finalamp, 'Phase stability', 'Amplitude stability', plotname) + + bestcycle = np.array(finalphase).argmin() + + self.writer.writerow(['phase', np.nan] + list(finalphase)) + self.writer.writerow(['amp', np.nan] + list(finalamp)) + + if len(finalphase) > 3: + phase_decrease, phase_quality, amp_quality = self.linreg_slope(finalphase[:4]), self.linreg_slope( + finalphase[-3:]), self.linreg_slope(finalamp[-3:]) + print(phase_decrease, phase_quality, amp_quality) + if phase_decrease < 0 and abs(phase_quality) < 0.05 and abs(amp_quality) < 0.05: + accept = True + else: + accept = False + return bestcycle, accept + else: + return None, None + + @staticmethod + def get_rms(fitsfile: str = None, maskSup: float = 1e-7): + """ + find the rms of an array, from Cycil Tasse/kMS + + :param fitsfile: fits file name + :param maskSup: mask threshold + + :return: rms --> rms of image + """ + hdul = fits.open(fitsfile) + mIn = np.ndarray.flatten(hdul[0].data) + m = mIn[np.abs(mIn) > maskSup] + rmsold = np.std(m) + diff = 1e-1 + cut = 3. + med = np.median(m) + for i in range(10): + ind = np.where(np.abs(m - med) < rmsold * cut)[0] + rms = np.std(m[ind]) + if np.abs((rms - rmsold) / rmsold) < diff: break + rmsold = rms + hdul.close() + + print(f'rms: {rms}') + + return rms # jy/beam + + @staticmethod + def get_minmax(fitsfile: str = None): + """ + Get min/max value + + :param fitsfile: fits file name + + :return: minmax --> pixel min/max value + """ + hdul = fits.open(fitsfile) + data = hdul[0].data + minmax = np.abs(data.min() / data.max()) + hdul.close() + print(f"min/max: {minmax}") + + return minmax + + def image_stability(self): + """ + Determine image stability + + :return: bestcycle --> best solution cycle + accept --> accept this selfcal + """ + + # cycles = [self.get_cycle_num(fts) for fts in self.fitsfiles] + rmss = [self.get_rms(fts) * 1000 for fts in self.fitsfiles] + minmaxs = [self.get_minmax(fts) for fts in self.fitsfiles] + entropy_images = [self.image_entropy(fts) for fts in self.fitsfiles] + entropy_models = [self.image_entropy(fts) for fts in self.modelfiles] + + self.make_figure(rmss, minmaxs, '$RMS (mJy)$', '$|min/max|$', f'image_stability_{self.sourcename}.png') + self.make_figure(vals1=entropy_images, vals2=entropy_models, label1='Entropy image', label2='Entropy model', + plotname=f'entropy_{self.sourcename}.png') + + self.writer.writerow(['min/max'] + minmaxs + [np.nan]) + self.writer.writerow(['rms'] + rmss + [np.nan]) + self.writer.writerow(['entropy_image'] + entropy_images + [np.nan]) + self.writer.writerow(['entropy_model'] + entropy_models + [np.nan]) + + # SCORING + best_rms_cycle, best_minmax_cycle, best_entropy_cycle = (np.array(rmss[1:]).argmin() + 1, + np.array(minmaxs[1:]).argmin() + 1, + min(np.array(entropy_images).argmin(), + np.array(entropy_models).argmin()) + 1) + # using maxmin instead of minmax due to easier slope value to work with + rms_slope, maxmin_slope = linregress(list(range(len(rmss))), rmss).slope, linregress(list(range(len(rmss))), + 1 / np.array( + minmaxs)).slope + + # accept direction or not + if maxmin_slope > 10: + accept = True + elif maxmin_slope < 0 and rms_slope > 0: + accept = False + elif maxmin_slope < 1.5: + accept = False + elif maxmin_slope < 2 and rms_slope >= 0: + accept = False + else: + accept = True + + # choose best cycle + if best_rms_cycle == best_minmax_cycle: + bestcycle = best_rms_cycle + elif rmss[best_minmax_cycle] <= 1.1: + bestcycle = best_minmax_cycle + elif minmaxs[best_rms_cycle] / minmaxs[0] <= 1: + bestcycle = best_rms_cycle + else: + bestcycle = int(round(np.mean([best_minmax_cycle, best_rms_cycle, best_entropy_cycle]), 0)) + + return bestcycle - 1, accept + + +def parse_args(): + """ + Command line argument parser + + :return: parsed arguments + """ + parser = ArgumentParser(description='Determine selfcal quality') + parser.add_argument('--selfcal_folder', default='.') + parser.add_argument('--remote_only', action='store_true', help='Only remote stations are considered', default=None) + parser.add_argument('--international_only', action='store_true', help='Only international stations are considered', + default=None) + return parser.parse_args() + + +def main(): + """ + Main function + """ + args = parse_args() + + sq = SelfcalQuality(args.selfcal_folder, args.remote_only, args.international_only) + bestcycle_solutions, accept_solutions = sq.solution_stability() + bestcycle_image, accept_image = sq.image_stability() + sq.textfile.close() + df = pd.read_csv('selfcal_performance.csv').set_index('solutions').T + print(df) + df.to_csv(f'selfcal_performance_{sq.sourcename}.csv', index=False) + + print(f"Best cycle according to solutions {bestcycle_solutions} (SCORES IN TESTING PHASE)") + print(f"Accept according to solutions {accept_solutions} (SCORES IN TESTING PHASE)") + print(f"Best cycle according to image {bestcycle_image} (SCORES IN TESTING PHASE)") + print(f"Accept according to image {accept_image} (SCORES IN TESTING PHASE)") + + +if __name__ == '__main__': + main() diff --git a/h5_helpers/smooth_bandpass.py b/h5_helpers/smooth_bandpass.py index 0c25b690..6031af6e 100644 --- a/h5_helpers/smooth_bandpass.py +++ b/h5_helpers/smooth_bandpass.py @@ -16,73 +16,101 @@ import os from argparse import ArgumentParser -tables.file._open_files.close_all() - -parser = ArgumentParser() -parser.add_argument('--antennas', nargs='+', help='Antenna indices', required=True) -parser.add_argument('--no_plot', action='store_true', help='Do not make plot at the end', default=None) -parser.add_argument('--h5', type=str, help='h5 table', required=True ) -args = parser.parse_args() - -change_antennas = [int(a) for a in args.antennas] -# change_antennas = [56, 57, 58, 73] # indices of antennas to fix def overwrite_val(T, table, values, title=None): + """ + Overwrite value table + """ + ss = T.root._f_get_child('calibrator') ss._f_get_child(table).val._f_remove() T.create_array(ss._f_get_child(table), 'val', values) + def flatlist(l): + """ + Flatten a list + """ + return [x for xs in l for x in xs] -newfile = args.h5.replace(".h5", "_smooth.h5") - -os.system(f'cp {args.h5} {newfile}') -H = tables.open_file(newfile, 'r+') - -freq = H.root.calibrator.bandpass.freq[:] - -new_bandpass = np.zeros(H.root.calibrator.bandpass.val.shape) - -N=0 -shape = H.root.calibrator.bandpass.val.shape -T=shape[0]*shape[2]*shape[3] -for bb_pol in range(len(H.root.calibrator.bandpass.pol)): - for time in range(len(H.root.calibrator.bandpass.time)): - for antenna in range(len(H.root.calibrator.bandpass.ant)): - bp = H.root.calibrator.bandpass.val[time, :, antenna, bb_pol] - if antenna in change_antennas: - bp_new = bp.copy() - outliers = np.abs(np.subtract(bp,medfilt(bp, 55)))/np.nanmean(medfilt(bp, 55))>30/np.nanmean(medfilt(bp, 55)) - outliers = sorted(set(flatlist([[j for j in list(range(m - 7, m + 7)) if j 30 / np.nanmean( + medfilt(bp, 55)) + outliers = sorted(set(flatlist([[j for j in list(range(m - 7, m + 7)) if j < len(outliers)] for m in + [i[0] for i in np.argwhere(outliers)]]))) + bp_new[outliers] = medfilt(bp, 55)[outliers] + new_bandpass[time, :, antenna, bb_pol] = bp_new + else: + new_bandpass[time, :, antenna, bb_pol] = bp + N += 1 + print(f'{int(100 * N / T)}/{100}%', end='\r') + + overwrite_val(H, 'bandpass', new_bandpass) + + if not args.no_plot: + T = tables.open_file(args.h5) + + fig = plt.figure(figsize=(10, 10)) + rows = len(change_antennas) // 2 + columns = 2 + grid = plt.GridSpec(rows, columns, wspace=.25, hspace=.25) + for i in range(rows * columns): + exec(f"plt.subplot(grid{[i]})") + plt.scatter(freq / 1000000, T.root.calibrator.bandpass.val[0, :, change_antennas[i], 0], s=1) + plt.scatter(freq / 1000000, H.root.calibrator.bandpass.val[0, :, change_antennas[i], 0], s=1) + plt.xlabel('freq [MHz]', fontsize=16) + plt.ylabel('amplitude', fontsize=16) + plt.xticks(size=14) + plt.yticks(size=14) + plt.legend(['old', 'new'], fontsize=14) + plt.title('Antenna index ' + str(change_antennas[i]), size=14) + plt.show() + + T.close() + + H.close() + + +if __name__ == "__main__": + main() diff --git a/h5_merger.py b/h5_merger.py index 0e30a692..289c0dfc 100644 --- a/h5_merger.py +++ b/h5_merger.py @@ -1,8 +1,13 @@ """ -USE THIS SCRIPT PREFERABLY WITH PYTHON 3, AS ALL TESTS ARE DONE WITH PYTHON 3. +h5parm merger script for merging h5parms containing phase and/or amplitude solutions for calibrating LOFAR observations. -For examples to run this script from the command line: +This script is a standalone script, such that it is easy to copy-paste it to other pipelines. + +For examples to run this script from the command line, please have a look at: https://github.com/jurjen93/lofar_helpers/blob/master/h5_merger_examples.md + +It is a work in progress, so for any bug reports or suggestions for improvements: +please contact jurjendejong@strw.leidenuniv.nl (or find me on LOFAR slack pages) """ __author__ = "Jurjen de Jong (jurjendejong@strw.leidenuniv.nl)" @@ -20,6 +25,7 @@ import sys import tables import warnings +from argparse import ArgumentParser warnings.filterwarnings('ignore') @@ -91,6 +97,7 @@ debug_message = 'Please contact jurjendejong@strw.leidenuniv.nl.' + def remove_numbers(inp): """ Remove numbers from string (keep only letters) @@ -100,6 +107,7 @@ def remove_numbers(inp): return "".join(re.findall("[a-zA-z]+", inp)) + def make_utf8(inp): """ Convert input to utf8 instead of bytes @@ -143,7 +151,7 @@ def overwrite_table(T, solset, table, values, title=None): title = 'Antenna names and positions' values = array(values, dtype=[('name', 'S16'), ('position', ' 0: print('Ignore MS for time and freq axis, as --h5_time_freq is given.') self.ax_time = array([]) @@ -295,18 +310,18 @@ def __init__(self, h5_out, h5_tables=None, ms_files=None, h5_time_freq=None, con else: print('No freq axes in ' + h5 + '/' + solset + '/' + soltab) h5.close() - elif type(h5_time_freq)==str: - if len(self.ms)>0: + elif type(h5_time_freq) == str: + if len(self.ms) > 0: print('Ignore MS for time and freq axis, as --h5_time_freq is given.') - print('Take the time and freq from the following h5 solution file:\n'+h5_time_freq) + print('Take the time and freq from the following h5 solution file:\n' + h5_time_freq) T = tables.open_file(h5_time_freq) self.ax_time = T.root.sol000.phase000.time[:] self.ax_freq = T.root.sol000.phase000.freq[:] T.close() # use ms files for available information - elif len(self.ms)>0: - print('Take the time and freq from the following measurement sets:\n'+'\n'.join(self.ms)) + elif len(self.ms) > 0: + print('Take the time and freq from the following measurement sets:\n' + '\n'.join(self.ms)) self.ax_time = array([]) self.ax_freq = array([]) for m in self.ms: @@ -322,7 +337,8 @@ def __init__(self, h5_out, h5_tables=None, ms_files=None, h5_time_freq=None, con # if no ms files, use the longest time and frequency resolution from h5 tables else: - print('No MS or h5 file given for time/freq axis.\nUse frequency and time axis by combining all input h5 tables.') + print( + 'No MS or h5 file given for time/freq axis.\nUse frequency and time axis by combining all input h5 tables.') self.ax_time = array([]) self.ax_freq = array([]) for h5_name in self.h5_tables: @@ -336,19 +352,18 @@ def __init__(self, h5_out, h5_tables=None, ms_files=None, h5_time_freq=None, con time = st._f_get_child('time')[:] self.ax_time = sort(unique(append(self.ax_time, time))) else: - print('No time axes in '+h5+'/'+solset+'/'+soltab) + print('No time axes in ' + h5 + '/' + solset + '/' + soltab) if 'freq' in axes: freq = st._f_get_child('freq')[:] self.ax_freq = sort(unique(append(self.ax_freq, freq))) else: - print('No freq axes in '+h5+'/'+solset+'/'+soltab) + print('No freq axes in ' + h5 + '/' + solset + '/' + soltab) h5.close() - # get polarization output axis and check number of error and tec tables in merge list self.polarizations, polarizations = [], [] self.doublefulljones = False - self.tecnum, self.errornum = 0, 0 # to average in case of multiple tables + self.tecnum, self.errornum = 0, 0 # to average in case of multiple tables for n, h5_name in enumerate(self.h5_tables): h5 = tables.open_file(h5_name) if 'phase000' in h5.root.sol000._v_children.keys() \ @@ -362,8 +377,8 @@ def __init__(self, h5_out, h5_tables=None, ms_files=None, h5_time_freq=None, con self.fulljones_amplitudes = OrderedDict() # take largest polarization list/array - if len(polarizations)>len(self.polarizations): - if type(self.polarizations)==list: + if len(polarizations) > len(self.polarizations): + if type(self.polarizations) == list: self.polarizations = polarizations else: self.polarizations = polarizations.copy() @@ -376,14 +391,14 @@ def __init__(self, h5_out, h5_tables=None, ms_files=None, h5_time_freq=None, con h5.close() # check if fulljones - if len(self.polarizations)==4: + if len(self.polarizations) == 4: self.fulljones = True - elif len(self.polarizations)>4 or len(self.polarizations)==3: + elif len(self.polarizations) > 4 or len(self.polarizations) == 3: sys.exit('Output solutions cannot have ' + str(len(self.polarizations)) + ' solutions') else: self.fulljones = False - print('Output polarization:\n'+str(self.polarizations)) + print('Output polarization:\n' + str(self.polarizations)) self.poldim = len(self.polarizations) # validation checks @@ -407,12 +422,11 @@ def __init__(self, h5_out, h5_tables=None, ms_files=None, h5_time_freq=None, con # directions in an ordered dictionary self.directions = OrderedDict() - if len(self.directions)>1 and self.doublefulljones: - sys.exit("ERROR: Merging not compatitable with multiple directions and double fuljones merge") #TODO: update + if len(self.directions) > 1 and self.doublefulljones: + sys.exit("ERROR: Merging not compatitable with multiple directions and double fuljones merge") # TODO: update self.merge_diff_freq = merge_diff_freq - @property def have_same_antennas(self): """ @@ -430,10 +444,10 @@ def have_same_antennas(self): for soltab1 in ss1._v_groups.keys(): if (len(antennas_ref['name']) != len(ss1._f_get_child(soltab1).ant[:])) or \ (not all(antennas_ref['name'] == ss1._f_get_child(soltab1).ant[:])): - message = '\n'.join(['\nMismatch in antenna tables in '+h5_name1, - 'Antennas from '+'/'.join([solset1, 'antenna']), + message = '\n'.join(['\nMismatch in antenna tables in ' + h5_name1, + 'Antennas from ' + '/'.join([solset1, 'antenna']), antennas_ref['name'], - 'Antennas from '+'/'.join([solset1, soltab1, 'ant']), + 'Antennas from ' + '/'.join([solset1, soltab1, 'ant']), ss1._f_get_child(soltab1).ant[:]]) print(message) H_ref.close() @@ -458,11 +472,12 @@ def have_same_antennas(self): antennas = ss2.antenna[:] if (len(antennas_ref['name']) != len(antennas['name'])) \ or (not all(antennas_ref['name'] == antennas['name'])): - message = '\n'.join(['\nMismatch between antenna tables from '+h5_name1+' and '+h5_name2, - 'Antennas from '+h5_name1+':', - antennas_ref['name'], - 'Antennas from '+h5_name2+':', - antennas['name']]) + message = '\n'.join( + ['\nMismatch between antenna tables from ' + h5_name1 + ' and ' + h5_name2, + 'Antennas from ' + h5_name1 + ':', + antennas_ref['name'], + 'Antennas from ' + h5_name2 + ':', + antennas['name']]) print(message) H.close() H_ref.close() @@ -514,7 +529,8 @@ def _unpack_h5(self, st, solset, soltab): print("No {av} in {solset}/{soltab}".format(av=av, solset=solset, soltab=soltab)) if len(st.getAxesNames()) != len(st.getValues()[0].shape): - sys.exit('ERROR: Axes ({axlen}) and Value dimensions ({vallen}) are not equal'.format(axlen=len(st.getAxesNames()), vallen=len(st.getValues()[0].shape))) + sys.exit('ERROR: Axes ({axlen}) and Value dimensions ({vallen}) are not equal'.format( + axlen=len(st.getAxesNames()), vallen=len(st.getValues()[0].shape))) print('Value shape before --> {values}'.format(values=st.getValues()[0].shape)) @@ -530,7 +546,8 @@ def _unpack_h5(self, st, solset, soltab): axes_current.insert(0, 'dir') else: axes_current.insert(1, 'dir') - values = reorderAxes(origin_values.reshape(origin_values.shape+(1,)), st.getAxesNames()+['dir'], axes_current) + values = reorderAxes(origin_values.reshape(origin_values.shape + (1,)), st.getAxesNames() + ['dir'], + axes_current) # current axes for reordering of axes self.axes_current = [an for an in self.solaxnames if an in st.getAxesNames()] @@ -544,7 +561,7 @@ def _unpack_h5(self, st, solset, soltab): # remove invalid values values = self.remove_invalid_values(st.getType(), values, self.axes_current) - if remove_numbers(st.getType())=='tec': + if remove_numbers(st.getType()) == 'tec': # add frequencies if 'freq' not in st.getAxesNames(): ax = self.axes_final.index('freq') - len(self.axes_final) @@ -590,7 +607,7 @@ def _unpack_h5(self, st, solset, soltab): # interpolate time axis values = self._interp_along_axis(values, time_axes, self.ax_time, - self.axes_current.index('time')) + self.axes_current.index('time')) # make sure that the interpolation is performed well when merging tables with different frequencies if self.merge_diff_freq: @@ -615,7 +632,7 @@ def _unpack_h5(self, st, solset, soltab): values_tmp[:, :, :, :, -1, ...] = 1 values_tmp[:, :, :, :, 0, ...] = 1 for idx_old, f_v in enumerate(freq_axes): - idx_new = list([round(i/1_000_000, 1) for i in self.ax_freq]).index(round(f_v/1_000_000, 1)) + idx_new = list([round(i / 1_000_000, 1) for i in self.ax_freq]).index(round(f_v / 1_000_000, 1)) if self.axes_current.index(ax) == 0: values_tmp[idx_new, ...] = values[idx_old, ...] @@ -638,8 +655,7 @@ def _unpack_h5(self, st, solset, soltab): # interpolate freq axis values = self._interp_along_axis(values, freq_axes, self.ax_freq, - self.axes_current.index('freq')) - + self.axes_current.index('freq')) return values @@ -721,16 +737,16 @@ def _make_template_values(self, soltab): num_dir = max(len(self.directions), 1) - if 'amplitude' in soltab and len(self.polarizations)>0: + if 'amplitude' in soltab and len(self.polarizations) > 0: self.amplitudes = ones( (max(len(self.polarizations), 1), num_dir, len(self.ant), len(self.ax_freq), len(self.ax_time))) - if len(self.polarizations)==4: + if len(self.polarizations) == 4: self.amplitudes[1, ...] = 0. self.amplitudes[2, ...] = 0. else: self.amplitudes = ones((num_dir, len(self.ant), len(self.ax_freq), len(self.ax_time))) - if 'phase' in soltab and len(self.polarizations)>0: + if 'phase' in soltab and len(self.polarizations) > 0: self.phases = zeros( (max(len(self.polarizations), 1), num_dir, len(self.ant), len(self.ax_freq), len(self.ax_time))) elif 'rotation' in soltab: @@ -738,13 +754,15 @@ def _make_template_values(self, soltab): else: self.phases = zeros((num_dir, len(self.ant), len(self.ax_freq), len(self.ax_time))) - if 'tec' in soltab and not self.convert_tec and len(self.polarizations)>0: - self.tec = zeros((max(len(self.polarizations), 1), num_dir, len(self.ant), len(self.ax_freq), len(self.ax_time))) + if 'tec' in soltab and not self.convert_tec and len(self.polarizations) > 0: + self.tec = zeros( + (max(len(self.polarizations), 1), num_dir, len(self.ant), len(self.ax_freq), len(self.ax_time))) elif 'tec' in soltab and not self.convert_tec: self.tec = zeros((num_dir, len(self.ant), len(self.ax_freq), len(self.ax_time))) - if 'error' in soltab and len(self.polarizations)>0: - self.error = zeros((max(len(self.polarizations), 1), num_dir, len(self.ant), len(self.ax_freq), len(self.ax_time))) + if 'error' in soltab and len(self.polarizations) > 0: + self.error = zeros( + (max(len(self.polarizations), 1), num_dir, len(self.ant), len(self.ax_freq), len(self.ax_time))) elif 'error' in soltab: self.error = zeros((num_dir, len(self.ant), len(self.ax_freq), len(self.ax_time))) @@ -784,7 +802,7 @@ def _interp_along_axis(x, interp_from, interp_to, axis): # axis with length 1 if len(interp_from) == 1: new_vals = x - for _ in range(len(interp_to)-1): + for _ in range(len(interp_to) - 1): new_vals = append(new_vals, x, axis=axis) else: interp_vals = interp1d(interp_from, x, axis=axis, kind='nearest', fill_value='extrapolate') @@ -820,7 +838,7 @@ def get_model_h5(self, solset, soltab): if not self.convert_tec or (self.convert_tec and 'tec' not in soltab): self._make_template_values(soltab) self.axes_final = [an for an in self.solaxnames if an in st.getAxesNames()] - if len(self.polarizations)>1 and 'pol' not in st.getAxesNames(): + if len(self.polarizations) > 1 and 'pol' not in st.getAxesNames(): self.axes_final.insert(0, 'pol') elif 'tec' in soltab and self.convert_tec: for st_group in self.all_soltabs: @@ -892,7 +910,7 @@ def remove_invalid_values(soltab, values, axlist): values[(~isfinite(values))] = 0. values[0, ...][(~isfinite(values[0, ...])) | (values[0, ...] == 0.)] = 1. values[-1, ...][(~isfinite(values[-1, ...])) | (values[-1, ...] == 0.)] = 1. - elif axlist.index('pol')-len(axlist)==-1: + elif axlist.index('pol') - len(axlist) == -1: values[(~isfinite(values))] = 0. values[..., 0][(~isfinite(values[..., 0])) | (values[..., 0] == 0.)] = 1. values[..., -1][(~isfinite(values[..., -1])) | (values[..., -1] == 0.)] = 1. @@ -903,7 +921,6 @@ def remove_invalid_values(soltab, values, axlist): return values - @staticmethod def _expand_poldim(values, dim_pol, type, haspol): """ @@ -917,7 +934,7 @@ def _expand_poldim(values, dim_pol, type, haspol): :return: input values with extra polarization axis """ - if dim_pol==3 or dim_pol>4: + if dim_pol == 3 or dim_pol > 4: sys.exit('ERROR: Invalid number of polarizations.' '\nOnly 1, 2, or 4 are allowed. ' 'This corresponds to XX, (XY, YX,) YY or RR, (RL, LR,) LL.') @@ -928,7 +945,7 @@ def _expand_poldim(values, dim_pol, type, haspol): if values.ndim < 5 and not haspol: if type == 'amplitude': values_new = ones((dim_pol,) + values.shape) - if dim_pol==4: + if dim_pol == 4: values_new[1, ...] = 0 values_new[2, ...] = 0 elif type in ['phase', 'error', 'tec']: @@ -936,7 +953,7 @@ def _expand_poldim(values, dim_pol, type, haspol): else: sys.exit('ERROR: Only type in [amplitude, phase] allowed [%s requested].' % str(type)) for i in range(dim_pol): - if dim_pol==4: + if dim_pol == 4: if i in [0, 3]: values_new[i, ...] = values else: @@ -945,7 +962,7 @@ def _expand_poldim(values, dim_pol, type, haspol): elif values.shape[0] in [1, 2] and dim_pol in [2, 4] and haspol: if type == 'amplitude': values_new = ones((dim_pol,) + values.shape[1:]) - if dim_pol==4: + if dim_pol == 4: values_new[1, ...] = 0 values_new[2, ...] = 0 elif type in ['phase', 'error', 'tec']: @@ -997,7 +1014,7 @@ def merge_tables(self, solset=None, soltab=None): # get values from h5 file table_values = self._unpack_h5(st, solset, soltab) - for dir_idx in range(num_dirs): # loop over all directions in h5 + for dir_idx in range(num_dirs): # loop over all directions in h5 if self.filtered_dir != None: if len(self.filtered_dir) > 0 and dir_idx not in self.filtered_dir: @@ -1013,9 +1030,11 @@ def merge_tables(self, solset=None, soltab=None): dirs = ss.getSou() if 'dir' in list(dirs.keys())[0].lower() and list(dirs.keys())[0][-1].isnumeric(): dirs = OrderedDict(sorted(dirs.items())) - elif len(dirs)>1 and (sys.version_info.major == 2 or (sys.version_info.major == 3 and sys.version_info.minor<6)): - print('WARNING: Order of source directions from h5 table might not be ordered. This is an old Python issue.' - '\nSuggest to switch to Python 3.6 or higher') + elif len(dirs) > 1 and ( + sys.version_info.major == 2 or (sys.version_info.major == 3 and sys.version_info.minor < 6)): + print( + 'WARNING: Order of source directions from h5 table might not be ordered. This is an old Python issue.' + '\nSuggest to switch to Python 3.6 or higher') # coordinate list source_coords = dirs[list(dirs.keys())[dir_idx]] @@ -1023,14 +1042,14 @@ def merge_tables(self, solset=None, soltab=None): if self.merge_all_in_one and self.n == 1: idx = 0 print('Merging direction {:f},{:f} with previous direction'.format(*source_coords)) - if abs(self.directions['Dir00'][0])>0 and abs(self.directions['Dir00'][1])>0: - self.add_direction({'Dir00': source_coords}) # 0.0 coordinate bug + if abs(self.directions['Dir00'][0]) > 0 and abs(self.directions['Dir00'][1]) > 0: + self.add_direction({'Dir00': source_coords}) # 0.0 coordinate bug print('Adding new direction {:f},{:f}'.format(*source_coords)) elif any([array_equal(source_coords, list(sv)) for sv in self.directions.values()]): # Direction already exists, add to the existing solutions. print('Direction {:f},{:f} already exists. Adding to this direction.'.format(*source_coords)) # Matching on 5 decimals rounding - idx = list([[round(l[0],5), round(l[1],5)] for l in self.directions.values()]).\ + idx = list([[round(l[0], 5), round(l[1], 5)] for l in self.directions.values()]). \ index([round(source_coords[0], 5), round(source_coords[1], 5)]) else: # new direction if abs(source_coords[0]) > 0 and abs(source_coords[1]) > 0: @@ -1049,7 +1068,7 @@ def merge_tables(self, solset=None, soltab=None): if self.n > shape[dir_index]: shape[dir_index] = 1 self.phases = append(self.phases, zeros(shape), - axis=dir_index) # add clean phase to merge with + axis=dir_index) # add clean phase to merge with elif st.getType() == 'amplitude': shape = list(self.amplitudes.shape) dir_index = self.amplitudes.ndim - 4 @@ -1058,7 +1077,7 @@ def merge_tables(self, solset=None, soltab=None): if self.n > shape[dir_index]: shape[dir_index] = 1 self.amplitudes = append(self.amplitudes, ones(shape), - axis=dir_index) # add clean gain to merge with + axis=dir_index) # add clean gain to merge with elif st.getType() == 'tec' and not self.convert_tec: shape = list(self.tec.shape) dir_index = self.tec.ndim - 4 @@ -1067,7 +1086,7 @@ def merge_tables(self, solset=None, soltab=None): if self.n > shape[dir_index]: shape[dir_index] = 1 self.tec = append(self.tec, zeros(shape), - axis=dir_index) # add clean phase to merge with + axis=dir_index) # add clean phase to merge with elif st.getType() == 'error': dir_index = self.error.ndim - 4 if dir_index < 0: @@ -1075,10 +1094,11 @@ def merge_tables(self, solset=None, soltab=None): if self.n > shape[dir_index]: shape[dir_index] = 1 self.error = append(self.error, zeros(shape), - axis=dir_index) # add clean phase to merge with + axis=dir_index) # add clean phase to merge with # Add the solution table axis --> tec, phase, rotation, amplitude, or error - if (st.getType() == 'tec' and self.convert_tec) or st.getType() == 'phase' or st.getType() == 'rotation': + if ( + st.getType() == 'tec' and self.convert_tec) or st.getType() == 'phase' or st.getType() == 'rotation': # add values if self.fulljones and not self.doublefulljones: @@ -1107,7 +1127,7 @@ def merge_tables(self, solset=None, soltab=None): self.phases[idx, ...] += values[...] print(diagdiag_math) - elif self.doublefulljones: # save table for in later double fulljones merge + elif self.doublefulljones: # save table for in later double fulljones merge self.fulljones_phases.update({h5_name: values}) else: @@ -1116,22 +1136,22 @@ def merge_tables(self, solset=None, soltab=None): elif st.getType() == 'tec' and not self.convert_tec: # add values - if len(self.polarizations)>0: + if len(self.polarizations) > 0: # average tec - self.tec[:, idx, ...] += values[:, ...]/self.tecnum + self.tec[:, idx, ...] += values[:, ...] / self.tecnum else: # average tec - self.tec[idx, ...] += values[...]/self.tecnum + self.tec[idx, ...] += values[...] / self.tecnum elif st.getType() == 'error': # add values - if len(self.polarizations)>0: + if len(self.polarizations) > 0: # average error - self.error[:, idx, ...] += values[:, ...]/self.errornum + self.error[:, idx, ...] += values[:, ...] / self.errornum else: # average error - self.error[idx, ...] += values[...]/self.errornum + self.error[idx, ...] += values[...] / self.errornum elif st.getType() == 'amplitude': @@ -1140,20 +1160,22 @@ def merge_tables(self, solset=None, soltab=None): if self.fulljones and not self.doublefulljones: if 'pol' in st.getAxesNames() and not fulljones_done: - if st.getAxisLen('pol')==4: + if st.getAxisLen('pol') == 4: print(self.amplitudes[:, 0, 0, 0, 0]) - self.amplitudes[1, idx, ...] = self.amplitudes[0, idx, ...] * values[1, ...] # Axx * Bxy - self.amplitudes[2, idx, ...] = self.amplitudes[-1, idx, ...] * values[2, ...] # Ayy * Byx + self.amplitudes[1, idx, ...] = self.amplitudes[0, idx, ...] * values[ + 1, ...] # Axx * Bxy + self.amplitudes[2, idx, ...] = self.amplitudes[-1, idx, ...] * values[ + 2, ...] # Ayy * Byx fulljones_done = True print(diagfull_math) elif fulljones_done: - self.amplitudes[1, idx, ...] *= values[-1, ...] # Axy * Byy - self.amplitudes[2, idx, ...] *= values[0, ...] # Ayx * Bxx + self.amplitudes[1, idx, ...] *= values[-1, ...] # Axy * Byy + self.amplitudes[2, idx, ...] *= values[0, ...] # Ayx * Bxx print(fulldiag_math) - self.amplitudes[0, idx, ...] *= values[0, ...] # Axx * Bxx - self.amplitudes[-1, idx, ...] *= values[-1, ...] # Ayy * Byy + self.amplitudes[0, idx, ...] *= values[0, ...] # Axx * Bxx + self.amplitudes[-1, idx, ...] *= values[-1, ...] # Ayy * Byy elif not self.doublefulljones: if 'pol' in self.axes_current: @@ -1162,7 +1184,7 @@ def merge_tables(self, solset=None, soltab=None): self.amplitudes[idx, ...] *= values[...] print(diagdiag_math) - elif self.doublefulljones: # save table for doublefulljones merge + elif self.doublefulljones: # save table for doublefulljones merge self.fulljones_amplitudes.update({h5_name: values}) else: @@ -1185,21 +1207,20 @@ def reshape_dir(arr): shape.insert(1, 1) return arr.reshape(shape) - print(doublefulljones_math) tables = [] for h5, t in self.fulljones_phases.items(): if h5 in self.fulljones_amplitudes.keys(): - tables.append(self.fulljones_amplitudes[h5].astype(complex128) * exp(t*1j)) + tables.append(self.fulljones_amplitudes[h5].astype(complex128) * exp(t * 1j)) else: - out = exp(t*1j) - out[1]=0. - out[2]=0. + out = exp(t * 1j) + out[1] = 0. + out[2] = 0. tables.append(out) pol_ax = self.axes_final.index('pol') - output_table = self.amplitudes.astype(complex128) * exp(self.phases*1j) + output_table = self.amplitudes.astype(complex128) * exp(self.phases * 1j) print("Calculating...") # Set template matrix @@ -1224,10 +1245,10 @@ def reshape_dir(arr): Byy = reshape_dir(take(input_table, indices=[3], axis=pol_ax)) # Necessary to make sure Axx, .. Ayy are not overwritten during calculations - nAxx = Axx*Bxx + Axy*Byx - nAxy = Axy*Byy + Axx*Bxy - nAyx = Ayx*Bxx + Ayy*Byx - nAyy = Ayy*Byy + Ayx*Bxy + nAxx = Axx * Bxx + Axy * Byx + nAxy = Axy * Byy + Axx * Bxy + nAyx = Ayx * Bxx + Ayy * Byx + nAyy = Ayy * Byy + Ayx * Bxy # Set final new values Axx = nAxx.copy() @@ -1238,7 +1259,7 @@ def reshape_dir(arr): print('Done\n') # Construct final output matrix - if pol_ax==0: + if pol_ax == 0: output_table[0, ...] = Axx output_table[1, ...] = Axy output_table[2, ...] = Ayx @@ -1290,19 +1311,19 @@ def reorder_directions(self): H = tables.open_file(self.h5name_out, 'r+') for solset in H.root._v_groups.keys(): ss = H.root._f_get_child(solset) - #Not needed if only 1 source + # Not needed if only 1 source if len(ss.source[:]) == 1: H.close() return self - #Problem when source table is empty + # Problem when source table is empty elif len(ss.source[:]) == 0: H.close() - sys.exit('ERROR: No sources in output table '+'/'.join([solset, 'source'])) - #No reordering needed + sys.exit('ERROR: No sources in output table ' + '/'.join([solset, 'source'])) + # No reordering needed elif all(ss.source[:]['name'] == ss._f_get_child(list(ss._v_groups.keys())[0]).dir[:]): H.close() return self - #Reordering needed + # Reordering needed else: sources = array(sort(ss.source[:]), dtype=[('name', 'S128'), ('dir', ' {values}'.format(values=weights.shape)) solsetout.makeSoltab('phase', axesNames=self.axes_final, axesVals=axes_vals, vals=self.phases, weights=weights) elif 'amplitude' in soltab: self.amplitudes = self.remove_invalid_values('amplitude', self.amplitudes, self.axes_final) - if shape_axes_vals!=self.amplitudes.shape: - self.amplitudes=reorderAxes(self.amplitudes, self.axes_current, self.axes_final) + if shape_axes_vals != self.amplitudes.shape: + self.amplitudes = reorderAxes(self.amplitudes, self.axes_current, self.axes_final) weights = ones(self.amplitudes.shape) print('Value shape after --> {values}'.format(values=weights.shape)) solsetout.makeSoltab('amplitude', axesNames=self.axes_final, axesVals=axes_vals, vals=self.amplitudes, weights=weights) elif 'tec' in soltab and not self.convert_tec: self.tec = self.remove_invalid_values('tec', self.tec, self.axes_final) - if shape_axes_vals!=self.tec.shape: - self.tec=reorderAxes(self.tec, self.axes_current, self.axes_final) + if shape_axes_vals != self.tec.shape: + self.tec = reorderAxes(self.tec, self.axes_current, self.axes_final) weights = ones(self.tec.shape) print('Value shape after --> {values}'.format(values=weights.shape)) solsetout.makeSoltab('tec', axesNames=self.axes_final, @@ -1416,8 +1437,8 @@ def create_new_dataset(self, solset, soltab): vals=self.tec, weights=weights) elif 'error' in soltab: self.error = self.remove_invalid_values('error', self.error, self.axes_final) - if shape_axes_vals!=self.error.shape: - self.error=reorderAxes(self.error, self.axes_current, self.axes_final) + if shape_axes_vals != self.error.shape: + self.error = reorderAxes(self.error, self.axes_current, self.axes_final) weights = ones(self.error.shape) print('Value shape after --> {values}'.format(values=weights.shape)) solsetout.makeSoltab('error', axesNames=self.axes_final, @@ -1441,7 +1462,7 @@ def add_empty_directions(self, add_directions=None): return self h5 = h5parm(self.h5name_out, readonly=True) - filetemp = self.h5name_out.replace('.h5','_temph5merger.h5') + filetemp = self.h5name_out.replace('.h5', '_temph5merger.h5') h5_temp = h5parm(filetemp, readonly=False) solset = h5.getSolset(self.solset) solsettemp = h5_temp.makeSolset(self.solset) @@ -1464,7 +1485,7 @@ def add_empty_directions(self, add_directions=None): dir_index = solutiontable.getAxesNames().index('dir') new_shape = list(values.shape) last_idx = new_shape[dir_index] - new_idx = last_idx+len(add_directions)-1 + new_idx = last_idx + len(add_directions) - 1 new_shape[dir_index] = new_idx if 'phase' in st: @@ -1486,16 +1507,17 @@ def add_empty_directions(self, add_directions=None): values_new[:, :, :, :, 0:last_idx, ...] = values weights = ones(values_new.shape) - solsettemp.makeSoltab(remove_numbers(st), axesNames=list(axes.keys()), axesVals=list(axes.values()), vals=values_new, - weights=weights) + solsettemp.makeSoltab(remove_numbers(st), axesNames=list(axes.keys()), axesVals=list(axes.values()), + vals=values_new, + weights=weights) - print('Default directions added for '+self.solset+'/'+st) - print('Shape change: '+str(values.shape)+' ---> '+str(values_new.shape)) + print('Default directions added for ' + self.solset + '/' + st) + print('Shape change: ' + str(values.shape) + ' ---> ' + str(values_new.shape)) h5.close() h5_temp.close() - os.system('rm '+self.h5name_out +' && mv '+filetemp+' '+self.h5name_out) + os.system('rm ' + self.h5name_out + ' && mv ' + filetemp + ' ' + self.h5name_out) return self @@ -1511,7 +1533,7 @@ def change_pol(self, single=False, nopol=False): T = tables.open_file(self.h5name_out, 'r+') # DP3 axes order - if self.poldim==5 or single: + if self.poldim == 5 or single: output_axes = ['time', 'freq', 'ant', 'dir', 'pol'] else: output_axes = ['time', 'ant', 'dir', 'freq'] @@ -1521,17 +1543,17 @@ def change_pol(self, single=False, nopol=False): for soltab in ss._v_groups.keys(): st = ss._f_get_child(soltab) - if self.poldim==0: - print('No polarization in '+solset+'/'+soltab) + if self.poldim == 0: + print('No polarization in ' + solset + '/' + soltab) if nopol: continue elif single: T.create_array(st, 'pol', array([b'I'], dtype='|S2')) - elif self.poldim==1: + elif self.poldim == 1: st.pol._f_remove() T.create_array(st, 'pol', array([b'I'], dtype='|S2')) - elif self.poldim>1 and (single or nopol): + elif self.poldim > 1 and (single or nopol): for axes in ['val', 'weight']: if not all(st._f_get_child(axes)[..., 0] == \ st._f_get_child(axes)[..., -1]): @@ -1540,13 +1562,12 @@ def change_pol(self, single=False, nopol=False): '\nERROR: No polarization reduction will be done.' '\nERROR: Do not use --no_pol or --single_pol.') - for axes in ['val', 'weight']: - if self.poldim>0 and single: + if self.poldim > 0 and single: newval = st._f_get_child(axes)[:, :, :, :, 0:1] - elif nopol and self.poldim>0: + elif nopol and self.poldim > 0: newval = st._f_get_child(axes)[:, :, :, :, 0] - elif self.poldim==0 and single: + elif self.poldim == 0 and single: newval = expand_dims(st._f_get_child(axes)[:], len(st._f_get_child(axes).shape)) else: continue @@ -1562,7 +1583,7 @@ def change_pol(self, single=False, nopol=False): atomtype = tables.Float64Atom() old_axes = [ax for ax in st._f_get_child(axes).attrs['AXES'].decode('utf8').split(',') + ['pol'] - if ax in output_axes] + if ax in output_axes] idx = [old_axes.index(ax) for ax in output_axes] newval = transpose(newval, idx) st._f_get_child(axes)._f_remove() @@ -1574,7 +1595,6 @@ def change_pol(self, single=False, nopol=False): print('Value shape after changing poldim--> ' + str(st._f_get_child('val')[:].shape)) - T.close() return self @@ -1584,7 +1604,7 @@ def add_h5_antennas(self): Add antennas to output table from H5 list. """ - print('Add antenna table from '+self.h5_tables[0]) + print('Add antenna table from ' + self.h5_tables[0]) T = tables.open_file(self.h5_tables[0]) antennas = T.root.sol000.antenna[:] T.close() @@ -1602,7 +1622,7 @@ def add_ms_antennas(self, keep_h5_interstations=None): :param keep_h5_interstations: keep h5 international stations """ - print('Add antenna table from '+self.ms[0]) + print('Add antenna table from ' + self.ms[0]) if len(self.ms) == 0: sys.exit("ERROR: Measurement set needed to add antennas. Use --ms.") @@ -1631,11 +1651,12 @@ def add_ms_antennas(self, keep_h5_interstations=None): antenna_index = attrsaxes.decode('utf8').split(',').index('ant') h5_antlist = [v.decode('utf8') for v in list(st.ant[:])] - if keep_h5_interstations: # keep international stations if these are not in MS + if keep_h5_interstations: # keep international stations if these are not in MS new_antlist = [station for station in ms_antlist if 'CS' in station] + \ [station for station in h5_antlist if 'ST' not in station] all_antennas = [a for a in unique(append(ms_antennas, h5_antennas), axis=0) if a[0] != 'ST001'] - antennas_new = [all_antennas[[a[0].decode('utf8') for a in all_antennas].index(na)] for na in new_antlist] # sorting + antennas_new = [all_antennas[[a[0].decode('utf8') for a in all_antennas].index(na)] for na in + new_antlist] # sorting else: new_antlist = ms_antlist antennas_new = ms_antennas @@ -1648,11 +1669,12 @@ def add_ms_antennas(self, keep_h5_interstations=None): superstation_index = h5_antlist.index('ST001') superstation = True except ValueError: - superstation=False + superstation = False print('No super station (ST001) in antennas.') for axes in ['val', 'weight']: - assert axes in list(st._v_children.keys()), axes+' not in .root.'+solset+'.'+soltab+' (not in axes)' + assert axes in list( + st._v_children.keys()), axes + ' not in .root.' + solset + '.' + soltab + ' (not in axes)' h5_values = st._f_get_child(axes)[:] shape = list(h5_values.shape) shape[antenna_index] = len(new_antlist) @@ -1671,7 +1693,7 @@ def add_ms_antennas(self, keep_h5_interstations=None): ms_values[:, :, :, idx, ...] += h5_values[:, :, :, idx_h5, ...] elif antenna_index == 4: ms_values[:, :, :, :, idx, ...] += h5_values[:, :, :, :, idx_h5, ...] - elif 'CS' in antenna and superstation: # core stations + elif 'CS' in antenna and superstation: # core stations if antenna_index == 0: ms_values[idx, ...] += h5_values[superstation_index, ...] elif antenna_index == 1: @@ -1682,9 +1704,9 @@ def add_ms_antennas(self, keep_h5_interstations=None): ms_values[:, :, :, idx, ...] += h5_values[:, :, :, superstation_index, ...] elif antenna_index == 4: ms_values[:, :, :, :, idx, ...] += h5_values[:, :, :, :, superstation_index, ...] - elif antenna not in h5_antlist and ('amplitude' in soltab or axes=='weight') \ + elif antenna not in h5_antlist and ('amplitude' in soltab or axes == 'weight') \ and 'RS' not in antenna and 'CS' not in antenna: - if axes=='val': + if axes == 'val': print('Add ' + antenna + ' to output H5 from MS') if antenna_index == 0: ms_values[idx, ...] = 1 @@ -1697,7 +1719,6 @@ def add_ms_antennas(self, keep_h5_interstations=None): elif antenna_index == 4: ms_values[:, :, :, :, idx, ...] = 1 - valtype = str(st._f_get_child(axes).dtype) if '16' in valtype: atomtype = tables.Float16Atom() @@ -1711,7 +1732,7 @@ def add_ms_antennas(self, keep_h5_interstations=None): st._f_get_child(axes)._f_remove() H.create_array(st, axes, ms_values.astype(valtype), atom=atomtype) st._f_get_child(axes).attrs['AXES'] = attrsaxes - print('Value shape after --> '+str(st.val.shape)) + print('Value shape after --> ' + str(st.val.shape)) H.close() @@ -1727,7 +1748,7 @@ def add_template(self): if 'amplitude000' not in soltabs: if 'phase000' in soltabs: self.amplitudes = ones(H.root.sol000.phase000.val.shape) - if len(self.polarizations)==4: + if len(self.polarizations) == 4: self.amplitudes[..., 1] = 0. self.amplitudes[..., 2] = 0. self.create_new_dataset('sol000', 'amplitude') @@ -1755,8 +1776,9 @@ def flag_stations(self): antennas_in = list(st.ant[:]) axes = make_utf8(weight.attrs['AXES']).split(',') if 'tec' in soltab: - soltab='phase'+soltab[-3:] - axes_output = make_utf8(H.root._f_get_child(solset)._f_get_child(soltab).weight.attrs['AXES']).split(',') + soltab = 'phase' + soltab[-3:] + axes_output = make_utf8( + H.root._f_get_child(solset)._f_get_child(soltab).weight.attrs['AXES']).split(',') ant_index = axes.index('ant') weight = reorderAxes(weight, axes, [i for i in axes_output if i in axes]) for a in range(weight.shape[ant_index]): @@ -1792,7 +1814,7 @@ def add_weights(self): for solset in H.root._v_groups.keys(): ss = H.root._f_get_child(solset) for n, soltab in enumerate(ss._v_groups.keys()): - print(soltab+', from:') + print(soltab + ', from:') st = ss._f_get_child(soltab) shape = st.val.shape weight_out = ones(shape) @@ -1810,12 +1832,12 @@ def add_weights(self): weight = reorderAxes(weight, axes, [a for a in axes_new if a in axes]) # important to do the following in case the input tables are not all different directions - m = min(weight_out.shape[axes.index('dir')]-1, m) + m = min(weight_out.shape[axes.index('dir')] - 1, m) newvals = self._interp_along_axis(weight, st2.time[:], st.time[:], axes_new.index('time')) newvals = self._interp_along_axis(newvals, st2.freq[:], st.freq[:], axes_new.index('freq')) - if weight.ndim != weight_out.ndim: # not the same value shape + if weight.ndim != weight_out.ndim: # not the same value shape if 'pol' not in axes and 'pol' in axes_new: newvals = expand_dims(newvals, axis=axes_new.index('pol')) if newvals.shape[-1] != weight_out.shape[-1]: @@ -1829,17 +1851,17 @@ def add_weights(self): + input_h5 + ': ' + str(axes) + '\n axes from ' + self.h5name_out + ': ' + str(axes_new) + '.\n' + debug_message) - elif set(axes) == set(axes_new) and 'pol' in axes: # same axes + elif set(axes) == set(axes_new) and 'pol' in axes: # same axes pol_index = axes_new.index('pol') - if weight_out.shape[pol_index] == newvals.shape[pol_index]: # same pol numbers + if weight_out.shape[pol_index] == newvals.shape[pol_index]: # same pol numbers print(weight_out.shape, m, axes, len(self.h5_tables)) weight_out[:, :, :, m, ...] *= newvals[:, :, :, 0, ...] - else: # not the same polarization axis + else: # not the same polarization axis if newvals.shape[pol_index] != newvals.shape[-1]: sys.exit('ERROR: Upsampling of weights bug due to polarization axis mismatch.\n' + debug_message) - if newvals.shape[pol_index] == 1: # new values have only 1 pol axis + if newvals.shape[pol_index] == 1: # new values have only 1 pol axis for i in range(weight_out.shape[pol_index]): weight_out[:, :, :, m, i] *= newvals[:, :, :, 0, 0] elif newvals.shape[pol_index] == 2 and weight_out.shape[pol_index] == 1: @@ -1853,30 +1875,30 @@ def add_weights(self): else: sys.exit('ERROR: Upsampling of weights bug due to unexpected polarization mismatch.\n' + debug_message) - elif set(axes) == set(axes_new) and 'pol' not in axes: # same axes but no pol + elif set(axes) == set(axes_new) and 'pol' not in axes: # same axes but no pol dirind = axes_new.index('dir') if weight_out.shape[dirind] != newvals.shape[dirind] and newvals.shape[dirind] == 1: - if len(weight_out.shape)==4: + if len(weight_out.shape) == 4: weight_out[:, :, m, :] *= newvals[:, :, 0, :] - elif len(weight_out.shape)==5: + elif len(weight_out.shape) == 5: weight_out[:, :, :, m, :] *= newvals[:, :, :, 0, :] elif weight_out.shape[dirind] != newvals.shape[dirind]: - sys.exit('ERROR: Upsampling of weights because same direction exists multiple times in input h5 (verify and/or remove --propagate_flags)') + sys.exit( + 'ERROR: Upsampling of weights because same direction exists multiple times in input h5 (verify and/or remove --propagate_flags)') else: weight_out *= newvals else: sys.exit('ERROR: Upsampling of weights bug due to unexpected missing axes.\n axes from ' - + input_h5 +': '+str(axes) +'\n axes from ' - + self.h5name_out +': '+str(axes_new)+'.\n' - + debug_message) + + input_h5 + ': ' + str(axes) + '\n axes from ' + + self.h5name_out + ': ' + str(axes_new) + '.\n' + + debug_message) T.close() st.weight[:] = weight_out H.close() - print('\n') return self @@ -1893,9 +1915,9 @@ def format_tables(self): st = ss._f_get_child(soltab) dirs = st._f_get_child('dir')[:] if len(sources[:]) > len(dirs): - difference = list(set(sources)-set(dirs)) + difference = list(set(sources) - set(dirs)) - newdir = list(st._f_get_child('dir')[:])+difference + newdir = list(st._f_get_child('dir')[:]) + difference st._f_get_child('dir')._f_remove() H.create_array(st, 'dir', array(newdir).astype('|S5')) @@ -1908,12 +1930,12 @@ def format_tables(self): shape = list(newval.shape) for _ in difference: shape[dir_ind] = 1 - if 'amplitude' in soltab or axes=='weight': + if 'amplitude' in soltab or axes == 'weight': newval = append(newval, ones(shape), - axis=dir_ind) + axis=dir_ind) else: newval = append(newval, zeros(shape), - axis=dir_ind) + axis=dir_ind) valtype = str(st._f_get_child(axes).dtype) if '16' in valtype: @@ -1934,6 +1956,7 @@ def format_tables(self): return self + def _create_h5_name(h5_name): """ Correct such that output has always .h5 as extension @@ -1947,6 +1970,7 @@ def _create_h5_name(h5_name): h5_name += '.h5' return h5_name + def _change_solset(h5, solset_in, solset_out, delete=True, overwrite=True): """ This function is to change the solset numbers. @@ -1966,11 +1990,12 @@ def _change_solset(h5, solset_in, solset_out, delete=True, overwrite=True): print('Succesfully copied ' + solset_in + ' to ' + solset_out) if delete: H.root._f_get_child(solset_in)._f_remove(recursive=True) - print('Removed ' + solset_in +' in output') + print('Removed ' + solset_in + ' in output') H.close() return + def output_check(h5): """ Validate output by checking and comparing tables. @@ -1984,68 +2009,70 @@ def output_check(h5): H = tables.open_file(h5) - #check number of solset + # check number of solset assert len(list(H.root._v_groups.keys())) == 1, \ - 'More than 1 solset in '+str(list(H.root._v_groups.keys()))+'. Only 1 is allowed for h5_merger.py.' + 'More than 1 solset in ' + str(list(H.root._v_groups.keys())) + '. Only 1 is allowed for h5_merger.py.' for solset in H.root._v_groups.keys(): - #check sol00.. name - assert 'sol' in solset, solset+' is a wrong solset name, should be sol***' + # check sol00.. name + assert 'sol' in solset, solset + ' is a wrong solset name, should be sol***' ss = H.root._f_get_child(solset) - #check antennas + # check antennas antennas = ss.antenna - assert antennas.attrs.FIELD_0_NAME == 'name', 'No name in '+'/'.join([solset,'antenna']) - assert antennas.attrs.FIELD_1_NAME == 'position', 'No coordinate in '+'/'.join([solset,'antenna']) + assert antennas.attrs.FIELD_0_NAME == 'name', 'No name in ' + '/'.join([solset, 'antenna']) + assert antennas.attrs.FIELD_1_NAME == 'position', 'No coordinate in ' + '/'.join([solset, 'antenna']) - #check sources + # check sources sources = ss.source - assert sources.attrs.FIELD_0_NAME == 'name', 'No name in '+'/'.join([solset,'source']) - assert sources.attrs.FIELD_1_NAME == 'dir', 'No coordinate in '+'/'.join([solset,'source']) + assert sources.attrs.FIELD_0_NAME == 'name', 'No name in ' + '/'.join([solset, 'source']) + assert sources.attrs.FIELD_1_NAME == 'dir', 'No coordinate in ' + '/'.join([solset, 'source']) for soltab in ss._v_groups.keys(): st = ss._f_get_child(soltab) assert st.val.shape == st.weight.shape, \ - 'weight '+str(st.weight.shape)+' and values '+str(st.val.shape)+' do not have same shape' + 'weight ' + str(st.weight.shape) + ' and values ' + str(st.val.shape) + ' do not have same shape' - #check if pol and/or dir are missing + # check if pol and/or dir are missing for pd in ['pol', 'dir']: assert not (st.val.ndim == 5 and pd not in list(st._v_children.keys())), \ - '/'.join([solset, soltab, pd])+' is missing' + '/'.join([solset, soltab, pd]) + ' is missing' - #check if freq, time, and ant arrays are missing + # check if freq, time, and ant arrays are missing for fta in ['freq', 'time', 'ant']: assert fta in list(st._v_children.keys()), \ - '/'.join([solset, soltab, fta])+' is missing' + '/'.join([solset, soltab, fta]) + ' is missing' - #check if val and weight have AXES + # check if val and weight have AXES for vw in ['val', 'weight']: assert 'AXES' in st._f_get_child(vw).attrs._f_list("user"), \ - 'AXES missing in '+'/'.join([solset, soltab, vw]) + 'AXES missing in ' + '/'.join([solset, soltab, vw]) - #check if dimensions of values match with length of arrays + # check if dimensions of values match with length of arrays for ax_index, ax in enumerate(st.val.attrs['AXES'].decode('utf8').split(',')): assert st.val.shape[ax_index] == len(st._f_get_child(ax)[:]), \ - ax+' length is not matching with dimension from val in ' + '/'.join([solset, soltab, ax]) + ax + ' length is not matching with dimension from val in ' + '/'.join([solset, soltab, ax]) - #check if ant and antennas have equal sizes + # check if ant and antennas have equal sizes if ax == 'ant': assert len(antennas[:]) == len(st._f_get_child(ax)[:]), \ - '/'.join([solset, 'antenna'])+' and '+'/'.join([solset, soltab, ax])+ ' do not have same length' + '/'.join([solset, 'antenna']) + ' and ' + '/'.join( + [solset, soltab, ax]) + ' do not have same length' # check if dir and sources have equal sizes if ax == 'dir': assert len(sources[:]) == len(st._f_get_child(ax)[:]), \ - '/'.join([solset, 'source'])+' and '+'/'.join([solset, soltab, ax])+ ' do not have same length' + '/'.join([solset, 'source']) + ' and ' + '/'.join( + [solset, soltab, ax]) + ' do not have same length' - #check if phase and amplitude have same shapes + # check if phase and amplitude have same shapes for soltab1 in ss._v_groups.keys(): if ('phase' in soltab or 'amplitude' in soltab) and ('phase' in soltab1 or 'amplitude' in soltab1): st1 = ss._f_get_child(soltab1) assert st.val.shape == st1.val.shape, \ - '/'.join([solset, soltab, 'val']) + ' shape: '+str(st.weight.shape)+\ - '/'.join([solset, soltab1, 'val']) + ' shape: '+str(st1.weight.shape) + '/'.join([solset, soltab, 'val']) + ' shape: ' + str(st.weight.shape) + \ + '/'.join([solset, soltab1, 'val']) + ' shape: ' + str(st1.weight.shape) H.close() @@ -2053,6 +2080,7 @@ def output_check(h5): return True + class PolChange: """ This Python class helps to convert polarization from linear to circular or vice versa. @@ -2132,7 +2160,6 @@ def circ2lin(G): YX = 1j * (G[..., -1] - G[..., 0]) YY = (G[..., 0] + G[..., -1]) - if G.shape[-1] == 4: XX += (G[..., 2] + G[..., 1]) XY += 1j * (G[..., 2] - G[..., 1]) @@ -2166,13 +2193,12 @@ def add_polarization(values, dim_pol): :return: input values with extra polarization axis """ - values_new = ones(values.shape+(dim_pol,)) + values_new = ones(values.shape + (dim_pol,)) for i in range(dim_pol): values_new[..., i] = values return values_new - def create_template(self, soltab): """ Make template of the gains with only ones @@ -2195,16 +2221,16 @@ def create_template(self, soltab): values = reorderAxes(solutiontable.getValues()[0], solutiontable.getAxesNames(), self.axes_names[0:-1]) - self.G = ones(values.shape+(2,)).astype(complex128) + self.G = ones(values.shape + (2,)).astype(complex128) except: - sys.exit('ERROR: Received '+str(solutiontable.getAxesNames())+ + sys.exit('ERROR: Received ' + str(solutiontable.getAxesNames()) + ', but expect at least [time, freq, ant, dir] or [time, freq, ant, dir, pol]') self.axes_vals = {'time': solutiontable.getAxisValues('time'), - 'freq': solutiontable.getAxisValues('freq'), - 'ant': solutiontable.getAxisValues('ant'), - 'dir': solutiontable.getAxisValues('dir'), - 'pol': ['XX', 'XY', 'YX', 'YY']} + 'freq': solutiontable.getAxisValues('freq'), + 'ant': solutiontable.getAxisValues('ant'), + 'dir': solutiontable.getAxisValues('dir'), + 'pol': ['XX', 'XY', 'YX', 'YY']} break print('Value shape {soltab} before --> {shape}'.format(soltab=soltab, shape=self.G.shape)) @@ -2236,7 +2262,8 @@ def add_tec(self, solutiontable): axes_vals_tec.update({'pol': ['XX', 'XY', 'YX', 'YY']}) axes_vals_tec = [v[1] for v in sorted(axes_vals_tec.items(), key=lambda pair: self.axes_names.index(pair[0]))] - self.solsetout.makeSoltab('tec', axesNames=tec_axes_names, axesVals=axes_vals_tec, vals=tec, weights=ones(tec.shape)) + self.solsetout.makeSoltab('tec', axesNames=tec_axes_names, axesVals=axes_vals_tec, vals=tec, + weights=ones(tec.shape)) return self @@ -2258,24 +2285,28 @@ def create_new_gain_table(self, lin2circ, circ2lin): print('{ss}/{st} from {h5}'.format(ss=ss, st=st, h5=self.h5in_name)) if 'phase' in st: if 'pol' in solutiontable.getAxesNames(): - values = reorderAxes(solutiontable.getValues()[0], solutiontable.getAxesNames(), self.axes_names) + values = reorderAxes(solutiontable.getValues()[0], solutiontable.getAxesNames(), + self.axes_names) self.G *= exp(values * 1j) else: - values = reorderAxes(solutiontable.getValues()[0], solutiontable.getAxesNames(), self.axes_names[0:-1]) + values = reorderAxes(solutiontable.getValues()[0], solutiontable.getAxesNames(), + self.axes_names[0:-1]) self.G *= exp(self.add_polarization(values, 2) * 1j) elif 'amplitude' in st: if 'pol' in solutiontable.getAxesNames(): - values = reorderAxes(solutiontable.getValues()[0], solutiontable.getAxesNames(), self.axes_names) + values = reorderAxes(solutiontable.getValues()[0], solutiontable.getAxesNames(), + self.axes_names) self.G *= values else: - values = reorderAxes(solutiontable.getValues()[0], solutiontable.getAxesNames(), self.axes_names[0:-1]) + values = reorderAxes(solutiontable.getValues()[0], solutiontable.getAxesNames(), + self.axes_names[0:-1]) self.G *= self.add_polarization(values, 2) elif 'tec' in st: self.add_tec(solutiontable) else: - print("WARNING: didn't include {st} in this h5_merger.py version yet".format(st=st)+ + print("WARNING: didn't include {st} in this h5_merger.py version yet".format(st=st) + "\nPlease create a ticket on github if this needs to be changed") if n == 0: @@ -2300,22 +2331,25 @@ def create_new_gain_table(self, lin2circ, circ2lin): phase = angle(G_new) amplitude = abs(G_new) - #upsample weights - if phase.shape!=weights.shape: + # upsample weights + if phase.shape != weights.shape: new_weights = ones(phase.shape) for s in range(new_weights.shape[-1]): - if len(weights.shape)==4: + if len(weights.shape) == 4: new_weights[..., s] *= weights[:] - elif len(weights.shape)==5: + elif len(weights.shape) == 5: new_weights[..., s] *= weights[..., 0] weights = new_weights - self.axes_vals = [v[1] for v in sorted(self.axes_vals.items(), key=lambda pair: self.axes_names.index(pair[0]))] + self.axes_vals = [v[1] for v in + sorted(self.axes_vals.items(), key=lambda pair: self.axes_names.index(pair[0]))] - self.solsetout.makeSoltab('phase', axesNames=self.axes_names, axesVals=self.axes_vals, vals=phase, weights=weights) + self.solsetout.makeSoltab('phase', axesNames=self.axes_names, axesVals=self.axes_vals, vals=phase, + weights=weights) print('Created new phase solutions') - self.solsetout.makeSoltab('amplitude', axesNames=self.axes_names, axesVals=self.axes_vals, vals=amplitude, weights=weights) + self.solsetout.makeSoltab('amplitude', axesNames=self.axes_names, axesVals=self.axes_vals, vals=amplitude, + weights=weights) print('Created new amplitude solutions') # convert the polarization names, such that it is clear if the h5 is in circular or linear polarization @@ -2325,12 +2359,12 @@ def create_new_gain_table(self, lin2circ, circ2lin): st = ss.getSoltab(soltab) if 'pol' in st.getAxesNames(): pols = st.getAxisValues('pol') - if len(pols)==2: + if len(pols) == 2: if lin2circ: st.setAxisValues('pol', ['RR', 'LL']) elif circ2lin: st.setAxisValues('pol', ['XX', 'YY']) - if len(pols)==4: + if len(pols) == 4: if lin2circ: st.setAxisValues('pol', ['RR', 'RL', 'LR', 'LL']) elif circ2lin: @@ -2367,57 +2401,62 @@ def h5_check(h5): :param h5: h5 solution file """ - print('%%%%%%%%%%%%%%%%%%%%%%%%%\nSOLUTION FILE CHECK START\n%%%%%%%%%%%%%%%%%%%%%%%%%\n\nSolution file name:\n'+h5) + print( + '%%%%%%%%%%%%%%%%%%%%%%%%%\nSOLUTION FILE CHECK START\n%%%%%%%%%%%%%%%%%%%%%%%%%\n\nSolution file name:\n' + h5) H = tables.open_file(h5) solsets = list(H.root._v_groups.keys()) - print('\nFollowing solution sets in '+h5+ ':\n'+'\n'.join(solsets)) + print('\nFollowing solution sets in ' + h5 + ':\n' + '\n'.join(solsets)) for solset in solsets: ss = H.root._f_get_child(solset) soltabs = list(ss._v_groups.keys()) - print('\nFollowing solution tables in '+solset+':\n' + '\n'.join(soltabs)) - print('\nFollowing stations in '+solset+':\n'+', '.join([make_utf8(a) for a in list(ss.antenna[:]['name'])])) - print('\nFollowing sources in ' + solset + ':\n' + '\n'.join([make_utf8(a['name'])+'-->'+str([a['dir'][0], a['dir'][1]]) for a in list(ss.source[:])])) + print('\nFollowing solution tables in ' + solset + ':\n' + '\n'.join(soltabs)) + print('\nFollowing stations in ' + solset + ':\n' + ', '.join( + [make_utf8(a) for a in list(ss.antenna[:]['name'])])) + print('\nFollowing sources in ' + solset + ':\n' + '\n'.join( + [make_utf8(a['name']) + '-->' + str([a['dir'][0], a['dir'][1]]) for a in list(ss.source[:])])) for soltab in soltabs: st = ss._f_get_child(soltab) for table in ['val', 'weight']: axes = make_utf8(st._f_get_child(table).attrs["AXES"]) - print('\n'+'/'.join([solset, soltab, table])+' axes:\n' + + print('\n' + '/'.join([solset, soltab, table]) + ' axes:\n' + axes) - print('/'.join([solset, soltab, table])+' shape:\n' + + print('/'.join([solset, soltab, table]) + ' shape:\n' + str(st._f_get_child(table).shape)) - if table=='weight': + if table == 'weight': weights = st._f_get_child(table)[:] element_sum = 1 for w in weights.shape: element_sum *= w - flagged = round(sum(weights==0.)/element_sum*100, 2) - print('/'.join([solset, soltab, table])+' flagged:\n'+str(flagged)+'%') + flagged = round(sum(weights == 0.) / element_sum * 100, 2) + print('/'.join([solset, soltab, table]) + ' flagged:\n' + str(flagged) + '%') if 'pol' in axes: - print('/'.join([solset, soltab, 'pol'])+':\n'+','.join([make_utf8(p) for p in list(st._f_get_child('pol')[:])])) + print('/'.join([solset, soltab, 'pol']) + ':\n' + ','.join( + [make_utf8(p) for p in list(st._f_get_child('pol')[:])])) if 'time' in axes: time = st._f_get_child('time')[:] # print('/'.join([solset, soltab, 'time']) + ' start:\n' + str(time[0])) # print('/'.join([solset, soltab, 'time']) + ' end:\n' + str(time[-1])) - if len(st._f_get_child('time')[:])>1: - print('/'.join([solset, soltab, 'time']) + ' time resolution:\n' + str(diff(st._f_get_child('time')[:])[0])) + if len(st._f_get_child('time')[:]) > 1: + print('/'.join([solset, soltab, 'time']) + ' time resolution:\n' + str( + diff(st._f_get_child('time')[:])[0])) if 'freq' in axes: freq = st._f_get_child('freq')[:] # print('/'.join([solset, soltab, 'freq']) + ' start:\n' + str(freq[0])) # print('/'.join([solset, soltab, 'freq']) + ' end:\n' + str(freq[-1])) - if len(st._f_get_child('freq')[:])>1: - print('/'.join([solset, soltab, 'freq']) + ' resolution:\n' + str(diff(st._f_get_child('freq')[:])[0])) + if len(st._f_get_child('freq')[:]) > 1: + print('/'.join([solset, soltab, 'freq']) + ' resolution:\n' + str( + diff(st._f_get_child('freq')[:])[0])) H.close() print('\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%\nSOLUTION FILE CHECK FINISHED\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n') return + def _checknan_input(h5): """ Check h5 on nan or 0 values :param h5: h5parm file - - :return: 'Done' (to print) """ H = tables.open_file(h5) @@ -2425,7 +2464,7 @@ def _checknan_input(h5): try: amp = H.root.sol000.amplitude000.val[:] print('Amplitude:') - print(amp[(~isfinite(amp)) | (amp==0.)]) + print(amp[(~isfinite(amp)) | (amp == 0.)]) del amp except: print('H.root.amplitude000 does not exist') @@ -2433,14 +2472,15 @@ def _checknan_input(h5): try: phase = H.root.sol000.phase000.val[:] print('Phase:') - print(phase[(~isfinite(phase)) | (phase==0.)]) + print(phase[(~isfinite(phase)) | (phase == 0.)]) del phase except: print('H.root.sol000.phase000 does not exist') H.close() - return 'Done' + return + def _degree_to_radian(d): """ @@ -2451,7 +2491,8 @@ def _degree_to_radian(d): :return: return value in radian """ - return pi*d/180 + return pi * d / 180 + def move_source_in_sourcetable(h5, overwrite=False, dir_idx=None, dra_degrees=0, ddec_degrees=0): """ @@ -2464,7 +2505,7 @@ def move_source_in_sourcetable(h5, overwrite=False, dir_idx=None, dra_degrees=0, """ if not overwrite: - os.system('cp '+h5+' '+h5.replace('.h5', '_upd.h5')) + os.system('cp ' + h5 + ' ' + h5.replace('.h5', '_upd.h5')) h5 = h5.replace('.h5', '_upd.h5') H = tables.open_file(h5, 'r+') sources = H.root.sol000.source[:] @@ -2475,6 +2516,7 @@ def move_source_in_sourcetable(h5, overwrite=False, dir_idx=None, dra_degrees=0, return + def check_freq_overlap(h5_tables): """ Verify if the frequency bands between the h5 tables overlap @@ -2485,20 +2527,25 @@ def check_freq_overlap(h5_tables): for h52 in h5_tables: F = tables.open_file(h52) try: - if h51!=h52 and F.root.sol000.phase000.freq[:].max() < H.root.sol000.phase000.freq[:].min(): - print(f"WARNING: frequency bands between {h51} and {h52} do not overlap, you might want to use " - f"--merge_diff_freq to merge over different frequency bands") + if h51 != h52 and F.root.sol000.phase000.freq[:].max() < H.root.sol000.phase000.freq[:].min(): + print(f"WARNING: frequency bands between {h51} and {h52} do not overlap, you might want to use " + f"--merge_diff_freq to merge over different frequency bands") except: try: - if h51 != h52 and F.root.sol000.amplitude000.freq[:].max() < H.root.sol000.amplitude000.freq[:].min(): + if h51 != h52 and F.root.sol000.amplitude000.freq[:].max() < H.root.sol000.amplitude000.freq[ + :].min(): print(f"WARNING: frequency bands between {h51} and {h52} do not overlap, you might want to use " f"--merge_diff_freq to merge over different frequency bands") except: pass + return + + def check_time_overlap(h5_tables): """ Verify if the time slots between the h5 tables overlap + :param h5_tables: h5parm input tables """ for h51 in h5_tables: @@ -2506,29 +2553,40 @@ def check_time_overlap(h5_tables): for h52 in h5_tables: F = tables.open_file(h52) try: - if h51!=h52 and F.root.sol000.phase000.time[:].max() < H.root.sol000.phase000.time[:].min(): - print(f"WARNING: time slots between {h51} and {h52} do not overlap, might result in interpolation issues") + if h51 != h52 and F.root.sol000.phase000.time[:].max() < H.root.sol000.phase000.time[:].min(): + print( + f"WARNING: time slots between {h51} and {h52} do not overlap, might result in interpolation issues") except: try: - if h51 != h52 and F.root.sol000.amplitude000.time[:].max() < H.root.sol000.amplitude000.time[:].min(): + if h51 != h52 and F.root.sol000.amplitude000.time[:].max() < H.root.sol000.amplitude000.time[ + :].min(): print( f"WARNING: time slots between {h51} and {h52} do not overlap, might result in interpolation issues") except: pass + return + + def running_mean(nparray, avgfactor): """ Running mean over numpy axis + :param nparray: numpy array :param avgfactor: averaging factor + + :return: running mean """ cs = cumsum(insert(nparray, 0, 0)) + return (cs[avgfactor:] - cs[:-avgfactor]) / float(avgfactor) + def merge_h5(h5_out=None, h5_tables=None, ms_files=None, h5_time_freq=None, convert_tec=True, merge_all_in_one=False, lin2circ=False, circ2lin=False, add_directions=None, single_pol=None, no_pol=None, use_solset='sol000', filtered_dir=None, add_cs=None, add_ms_stations=None, check_output=None, freq_av=None, time_av=None, - check_flagged_station=True, propagate_flags=None, merge_diff_freq=None, no_antenna_crash=None, output_summary=None): + check_flagged_station=True, propagate_flags=None, merge_diff_freq=None, no_antenna_crash=None, + output_summary=None): """ Main function that uses the class MergeH5 to merge h5 tables. @@ -2558,7 +2616,7 @@ def merge_h5(h5_out=None, h5_tables=None, ms_files=None, h5_time_freq=None, conv tables.file._open_files.close_all() - if type(h5_tables)==str: + if type(h5_tables) == str: h5_tables = glob(h5_tables) if h5_out is None: @@ -2567,7 +2625,7 @@ def merge_h5(h5_out=None, h5_tables=None, ms_files=None, h5_time_freq=None, conv return print('\n##################################\nSTART MERGE HDF5 TABLES FOR LOFAR\n##################################' - '\n\nMerging the following tables:\n'+'\n'.join(h5_tables)+'\n') + '\n\nMerging the following tables:\n' + '\n'.join(h5_tables) + '\n') # Make sure that the h5 output file is not in the input list (as that will cause conflicts) if h5_out in h5_tables: @@ -2580,8 +2638,8 @@ def merge_h5(h5_out=None, h5_tables=None, ms_files=None, h5_time_freq=None, conv if use_solset != 'sol000': for h5_ind, h5 in enumerate(h5_tables): temph5 = h5.replace('.h5', '_temph5merger.h5') - print('Using different solset. Make temporary h5 file: '+temph5) - os.system('cp '+h5+' '+temph5) + print('Using different solset. Make temporary h5 file: ' + temph5) + os.system('cp ' + h5 + ' ' + temph5) _change_solset(temph5, use_solset, 'sol000') h5_tables[h5_ind] = temph5 @@ -2601,16 +2659,16 @@ def merge_h5(h5_out=None, h5_tables=None, ms_files=None, h5_time_freq=None, conv # Time averaging if time_av: - print("Time averaging with factor "+str(time_av)) - if len(merge.ax_time)>=time_av: + print("Time averaging with factor " + str(time_av)) + if len(merge.ax_time) >= time_av: merge.ax_time = running_mean(merge.ax_time, time_av)[::int(time_av)] else: print("WARNING: Time averaging factor larger than time axis, no averaging applied") # Freq averaging if freq_av: - print("Frequency averaging with factor "+str(freq_av)) - if len(merge.ax_freq)>=freq_av: + print("Frequency averaging with factor " + str(freq_av)) + if len(merge.ax_freq) >= freq_av: merge.ax_freq = running_mean(merge.ax_freq, freq_av)[::int(freq_av)] else: print("WARNING: Frequency averaging factor larger than frequency axis, no averaging applied") @@ -2638,13 +2696,12 @@ def merge_h5(h5_out=None, h5_tables=None, ms_files=None, h5_time_freq=None, conv # Make sure direction tables are in the same format merge.format_tables() - # Propagate weight flags from input into output if propagate_flags: merge.add_weights() # Add antennas - if (add_cs or add_ms_stations) and len(merge.ms)==0: + if (add_cs or add_ms_stations) and len(merge.ms) == 0: sys.exit('ERROR: --add_cs and --add_ms_stations need an MS, given with --ms.') if add_cs: merge.add_ms_antennas(keep_h5_interstations=True) @@ -2678,11 +2735,10 @@ def merge_h5(h5_out=None, h5_tables=None, ms_files=None, h5_time_freq=None, conv if no_pol: print('Remove polarization') merge.change_pol(nopol=True) - elif single_pol or merge.poldim==0: + elif single_pol or merge.poldim == 0: print('Make a single polarization') merge.change_pol(single=True) - ################################################# ############ POLARIZATION CONVERSION ############ ################################################# @@ -2693,10 +2749,10 @@ def merge_h5(h5_out=None, h5_tables=None, ms_files=None, h5_time_freq=None, conv elif lin2circ or circ2lin: if lin2circ: - h5_polchange = h5_out[0:-3]+'_circ.h5' + h5_polchange = h5_out[0:-3] + '_circ.h5' print('\nPolarization will be converted from linear to circular') else: - h5_polchange = h5_out[0:-3]+'_lin.h5' + h5_polchange = h5_out[0:-3] + '_lin.h5' print('\nPolarization will be converted from circular to linear') # Polarization conversion @@ -2709,18 +2765,18 @@ def merge_h5(h5_out=None, h5_tables=None, ms_files=None, h5_time_freq=None, conv Pol.create_new_gain_table(lin2circ, circ2lin) Pol.add_antenna_source_tables() - os.system('rm '+h5_out+' && cp '+h5_polchange+' '+h5_out+' && rm '+h5_polchange) + os.system('rm ' + h5_out + ' && cp ' + h5_polchange + ' ' + h5_out + ' && rm ' + h5_polchange) # Cleanup temp files for h5 in h5_tables: if '_temph5merger.h5' in h5: - os.system('rm '+h5) + os.system('rm ' + h5) # Validate output if check_output: output_check(h5_out) - print('\nSee output file --> '+h5_out+'\n\n###################\nEND MERGE H5 TABLES \n###################\n') + print('\nSee output file --> ' + h5_out + '\n\n###################\nEND MERGE H5 TABLES \n###################\n') try: # Close all tables that are still open @@ -2732,19 +2788,22 @@ def merge_h5(h5_out=None, h5_tables=None, ms_files=None, h5_time_freq=None, conv if output_summary: h5_check(h5_out) + return -if __name__ == '__main__': - from argparse import ArgumentParser - if sys.version_info.major == 2: - print('WARNING: This code is optimized for Python 3. Please switch to Python 3 if possible.') +def parse_input(): + """ + Command line argument parser - parser = ArgumentParser() + :return: arguments in correct format + """ + + parser = ArgumentParser(description='h5parm merger script for merging h5parms containing phase and/or amplitude solutions for calibrating LOFAR observations') parser.add_argument('-out', '--h5_out', type=str, help='Solution file output name. This name cannot be in the list of input solution files.') parser.add_argument('-in', '--h5_tables', type=str, nargs='+', help='Solution files to merge (can be both given as list with wildcard or string).', required=True) parser.add_argument('-ms', '--ms', type=str, help='Measurement Set input files (can be both given as list with wildcard or string).') parser.add_argument('--h5_time_freq', type=str, help='Solution file to use time and frequency arrays from. This is useful if the input solution files do not have the preferred time/freq resolution. ' - 'Input can be boolean (true/fals or 1/0) to use all h5 files or input can be 1 specific h5 file ().') + 'Input can be boolean (true/fals or 1/0) to use all h5 files or input can be 1 specific h5 file ().') parser.add_argument('--time_av', type=int, help='Time averaging factor.') parser.add_argument('--freq_av', type=int, help='Frequency averaging factor.') parser.add_argument('--keep_tec', action='store_true', help='Do not convert TEC to phase.') @@ -2766,12 +2825,11 @@ def merge_h5(h5_out=None, h5_tables=None, ms_files=None, h5_time_freq=None, conv parser.add_argument('--check_output', action='store_true', default=None, help='Check if the output has all the correct output information.') parser.add_argument('--merge_diff_freq', action='store_true', default=None, help='Merging tables over different frequency bands') - # parser.add_argument('--keep_sourcenames', action='store_true', default=None, help='Keep the name of the input sources') args = parser.parse_args() # check if solset name is accepted if 'sol' not in args.usesolset or sum([c.isdigit() for c in args.usesolset]) != 3: - sys.exit(args.usesolse+' not an accepted name. Only sol000, sol001, sol002, ... are accepted names for solsets.') + sys.exit(args.usesolse + ' not an accepted name. Only sol000, sol001, sol002, ... are accepted names for solsets.') if args.filter_directions: if (args.filter_directions.startswith("[") and args.filter_directions.endswith("]")): @@ -2781,59 +2839,67 @@ def merge_h5(h5_out=None, h5_tables=None, ms_files=None, h5_time_freq=None, conv sys.exit('--filter_directions can only have integers in the list.') else: filtered_dir[n] = int(v) + args.filter_directions = filtered_dir else: sys.exit('--filter_directions given but no list format. Please pass a list to --filter_directions.') else: - filtered_dir = [] + args.filter_directions = [] # make sure h5 tables in right format - if '[' in args.h5_tables: - h5tables = args.h5_tables.replace('[', '').replace(']', '').replace(' ', '').split(',') + if '[' in args.h5_tables[0]: + args.h5_tables = args.h5_tables[0].replace('[', '').replace(']', '').replace(' ', '').split(',') elif ' ' in args.h5_tables: - h5tables = args.h5_tables.split() - else: - h5tables = args.h5_tables + args.h5_tables = args.h5_tables[0].split() - if type(h5tables) == str: - h5tables = glob(h5tables) - elif type(h5tables) == list and len(h5tables) == 1: - h5tables = glob(h5tables[0]) - elif type(h5tables) == list: + if type(args.h5_tables) == str: + args.h5_tables = glob(args.h5_tables) + elif type(args.h5_tables) == list and len(args.h5_tables) == 1: + args.h5_tables = glob(args.h5_tables[0]) + elif type(args.h5_tables) == list: h5tablestemp = [] - for h5 in h5tables: + for h5 in args.h5_tables: h5tablestemp += glob(h5) + args.h5_tables = h5tablestemp if args.add_direction: - add_direction = args.add_direction.replace('[','').replace(']','').split(',') - add_direction = [float(add_direction[0]), float(add_direction[1])] - if add_direction[0] > pi*6 or add_direction[1] > pi*6: + args.add_direction = args.add_direction.replace('[', '').replace(']', '').split(',') + args.add_direction = [float(args.add_direction[0]), float(args.add_direction[1])] + if args.add_direction[0] > pi * 6 or args.add_direction[1] > pi * 6: sys.exit('ERROR: Please give --add_direction values in radian.') - else: - add_direction = None if args.h5_time_freq is not None: - if args.h5_time_freq.lower()=='true' or str(args.h5_time_freq)=='1': - h5_time_freq = True - elif args.h5_time_freq.lower()=='false' or str(args.h5_time_freq)=='0': - h5_time_freq = False - else: - h5_time_freq = args.h5_time_freq + if args.h5_time_freq.lower() == 'true' or str(args.h5_time_freq) == '1': + args.h5_time_freq = True + elif args.h5_time_freq.lower() == 'false' or str(args.h5_time_freq) == '0': + args.h5_time_freq = False else: - h5_time_freq = False + args.h5_time_freq = False + + return args + + +def main(): + """ + Main function + """ + if sys.version_info.major == 2: + print('WARNING: This code is optimized for Python 3. Please switch to Python 3 if possible.') + + args = parse_input() merge_h5(h5_out=args.h5_out, - h5_tables=h5tables, + h5_tables=args.h5_tables, ms_files=args.ms, - h5_time_freq=h5_time_freq, + h5_time_freq=args.h5_time_freq, convert_tec=not args.keep_tec, merge_all_in_one=args.merge_all_in_one, lin2circ=args.lin2circ, circ2lin=args.circ2lin, - add_directions=add_direction, + add_directions=args.add_direction, single_pol=args.single_pol, no_pol=args.no_pol, use_solset=args.usesolset, - filtered_dir=filtered_dir, + filtered_dir=args.filter_directions, add_cs=args.add_cs, add_ms_stations=args.add_ms_stations, check_output=args.check_output, @@ -2844,3 +2910,7 @@ def merge_h5(h5_out=None, h5_tables=None, ms_files=None, h5_time_freq=None, conv merge_diff_freq=args.merge_diff_freq, no_antenna_crash=args.no_antenna_crash, output_summary=args.output_summary) + + +if __name__ == '__main__': + main() diff --git a/make_boxes.py b/make_boxes.py index b9750e39..d19f25e6 100644 --- a/make_boxes.py +++ b/make_boxes.py @@ -36,64 +36,7 @@ warnings.filterwarnings("ignore") -parser = ArgumentParser() -parser.add_argument('-f', '--file', type=str, help='fitsfile name') -parser.add_argument('-l', '--location', type=str, help='data location folder name', default='.') -parser.add_argument('--no_images', action='store_true', help='store images') -parser.add_argument('--make_DL_food', action='store_true', help='store images for creating the DL model') -parser.add_argument('--ds9', action='store_true', help='open ds9 to interactively check and change the box selection') -parser.add_argument('-ac', '--angular_cutoff', type=float, default=None, help='angular distances higher than this value from the center will be excluded from the box selection') -parser.add_argument('-mb', '--max_boxes', type=int, default=999, help='Set max number of boxes that can be made') -args = parser.parse_args() -print(args) - -if sys.version_info.major == 2: - print('ERROR: This code only works for Python 3. Please switch.\n.....ENDED.....') - sys.exit() - -folder = args.location -if folder[-1] == '/': - folder = folder[0:-1] - -if args.make_DL_food: - os.system('mkdir -p {LOCATION}/DL_data/numpy'.format(LOCATION=folder)) - os.system('mkdir -p {LOCATION}/DL_data/png'.format(LOCATION=folder)) - if not os.path.isfile("{LOCATION}/DL_data/label_data.csv".format(LOCATION=folder)): - with open("{LOCATION}/DL_data/label_data.csv".format(LOCATION=folder), 'w') as file: - writer = csv.writer(file) - writer.writerow(["Name", "Recalibrate"]) - # try: # Install tkinter/tk for using interactive window mode - # import importlib - # importlib_found = importlib.util.find_spec("tk") is not None - # if not importlib_found: - # os.system('pip install --user tk') - # try: - # import tk - # use_tk = True - # except ModuleNotFoundError: - # use_tk = False - # except: - # print('ERROR: tk cannot be installed') - # use_tk = False - - -#check if folder exists and create if not -folders = folder.split('/') -for i, f in enumerate(folders): - subpath = '/'.join(folder.split('/')[0:i+1]) - if not os.path.isdir(subpath): - print(f'Create directory: {subpath}') - os.system(f'mkdir {subpath}') - -if not args.no_images: - try: - import matplotlib - matplotlib.use('Agg') - import matplotlib.pyplot as plt - from matplotlib.colors import SymLogNorm - except ImportError: - print('Failed to import matplotlib. Check your version.\nNo images will be made.') - args.no_images = True + def resample_pixels(image_data, rows, cols): """Resample image by summing pixels together""" @@ -495,8 +438,77 @@ def number_of_sources(self, sources): print(f'There are multiple sources in same box.\n Splitting possible.') return self +def parse_args(): + """Command line argument parser""" -if __name__ == '__main__': + parser = ArgumentParser() + parser.add_argument('-f', '--file', type=str, help='fitsfile name') + parser.add_argument('-l', '--location', type=str, help='data location folder name', default='.') + parser.add_argument('--no_images', action='store_true', help='store images') + parser.add_argument('--make_DL_food', action='store_true', help='store images for creating the DL model') + parser.add_argument('--ds9', action='store_true', + help='open ds9 to interactively check and change the box selection') + parser.add_argument('-ac', '--angular_cutoff', type=float, default=None, + help='angular distances higher than this value from the center will be excluded from the box selection') + parser.add_argument('-mb', '--max_boxes', type=int, default=999, help='Set max number of boxes that can be made') + args = parser.parse_args() + print(args) + + return args + +def main(): + """Main function""" + + args = parse_args() + + + if sys.version_info.major == 2: + print('ERROR: This code only works for Python 3. Please switch.\n.....ENDED.....') + sys.exit() + + folder = args.location + if folder[-1] == '/': + folder = folder[0:-1] + + if args.make_DL_food: + os.system('mkdir -p {LOCATION}/DL_data/numpy'.format(LOCATION=folder)) + os.system('mkdir -p {LOCATION}/DL_data/png'.format(LOCATION=folder)) + if not os.path.isfile("{LOCATION}/DL_data/label_data.csv".format(LOCATION=folder)): + with open("{LOCATION}/DL_data/label_data.csv".format(LOCATION=folder), 'w') as file: + writer = csv.writer(file) + writer.writerow(["Name", "Recalibrate"]) + # try: # Install tkinter/tk for using interactive window mode + # import importlib + # importlib_found = importlib.util.find_spec("tk") is not None + # if not importlib_found: + # os.system('pip install --user tk') + # try: + # import tk + # use_tk = True + # except ModuleNotFoundError: + # use_tk = False + # except: + # print('ERROR: tk cannot be installed') + # use_tk = False + + # check if folder exists and create if not + folders = folder.split('/') + for i, f in enumerate(folders): + subpath = '/'.join(folder.split('/')[0:i + 1]) + if not os.path.isdir(subpath): + print(f'Create directory: {subpath}') + os.system(f'mkdir {subpath}') + + if not args.no_images: + try: + import matplotlib + + matplotlib.use('Agg') + import matplotlib.pyplot as plt + from matplotlib.colors import SymLogNorm + except ImportError: + print('Failed to import matplotlib. Check your version.\nNo images will be made.') + args.no_images = True image = SetBoxes(fits_file=args.file) @@ -620,3 +632,6 @@ def number_of_sources(self, sources): if args.ds9: from box_helpers.move_boxes import move_boxes move_boxes(args.file, folder) + +if __name__ == '__main__': + main() diff --git a/ms_helpers/ms_flagger.py b/ms_helpers/ms_flagger.py index 52e68956..757dde97 100644 --- a/ms_helpers/ms_flagger.py +++ b/ms_helpers/ms_flagger.py @@ -1,13 +1,24 @@ from argparse import ArgumentParser import os -if __name__ == '__main__': +def parse_args(): + """ + Command line argument parser + :return: parsed arguments + """ parser = ArgumentParser() parser.add_argument('--ms', nargs='+', help='Input ms', required=True) parser.add_argument('--freqrange', help='Frequency range in MHz (example: 140-158)', required=True) parser.add_argument('--ant', help='Antenna (example: RS409HBA)', required=True) - args = parser.parse_args() + return parser.parse_args() + +def main(): + """ + Main function + """ + + args = parse_args() min_freq, max_freq = [float(f) for f in args.freqrange.split('-')] @@ -22,4 +33,7 @@ f"flag.freqrange='[{min_freq} .. {max_freq} MHz]'"] print(' '.join(command)) - os.system(' '.join(command)) \ No newline at end of file + os.system(' '.join(command)) + +if __name__ == '__main__': + main() diff --git a/phasediff_scores/phasediff_selection/phasediff_output.py b/phasediff_scores/phasediff_selection/phasediff_output.py index c133b594..00354699 100644 --- a/phasediff_scores/phasediff_selection/phasediff_output.py +++ b/phasediff_scores/phasediff_selection/phasediff_output.py @@ -8,6 +8,7 @@ from glob import glob import csv import sys +from argparse import ArgumentParser # sys.path.append("..") # from find_solint import GetSolint @@ -209,14 +210,25 @@ def theoretical_curve(self, t): self.C = self._get_C return limit * np.sqrt(1 - np.exp(-(self.C / (2 * t)))) -if __name__ == '__main__': - import argparse +def parse_args(): + """ + Command line argument parser - parser = argparse.ArgumentParser() + :return: parsed arguments + """ + parser = ArgumentParser() parser.add_argument('--h5', nargs='+', help='selfcal phasediff solutions', default=None) parser.add_argument('--station', help='for specific station', default=None) parser.add_argument('--all_stations', action='store_true', help='for all stations specifically') - args = parser.parse_args() + return parser.parse_args() + + + +if __name__ == '__main__': + + ###STILL AN EXPERIMENT!### + + args = parse_args() # set std score, for which you want to find the solint optimal_score = 2.5 diff --git a/subtract/subtract_with_dp3.py b/subtract/subtract_with_dp3.py index a9058cf7..53b117c9 100644 --- a/subtract/subtract_with_dp3.py +++ b/subtract/subtract_with_dp3.py @@ -319,10 +319,7 @@ def run(self, type=''): return self - - -if __name__ == "__main__": - +def main(): from argparse import ArgumentParser parser = ArgumentParser(description='Subtract region with WSClean') @@ -418,3 +415,6 @@ def run(self, type=''): concat=args.concat, applybeam=args.applybeam, applycal_h5=applycalh5, dirname=dirname) if not args.print_only_commands: Subtract.run(type='phaseshift') + +if __name__ == "__main__": + main() diff --git a/subtract/subtract_with_wsclean.py b/subtract/subtract_with_wsclean.py index 62c942f1..3ff15177 100644 --- a/subtract/subtract_with_wsclean.py +++ b/subtract/subtract_with_wsclean.py @@ -10,6 +10,7 @@ from itertools import repeat import re import pandas as pd +from argparse import ArgumentParser def add_trailing_zeros(s, digitsize=4): @@ -530,10 +531,12 @@ def run_DP3(self, phaseshift: str = None, freqavg: str = None, return self +def main(): + """ + Main function -if __name__ == "__main__": - - from argparse import ArgumentParser + :param args: Arguments parsed by ArgumentParser + """ parser = ArgumentParser(description='Subtract region with WSClean') parser.add_argument('--mslist', nargs='+', help='measurement sets', required=True) @@ -560,6 +563,7 @@ def run_DP3(self, phaseshift: str = None, freqavg: str = None, help='will search for the polygon_info.csv file to extract information from') parser.add_argument('--skip_predict', action='store_true', help='skip predict and do only subtract') parser.add_argument('--even_time_avg', action='store_true', help='(only if --forwidefield) only allow even time averaging (in case of stacking nights with different averaging)') + args = parser.parse_args() if not args.skip_predict: @@ -690,4 +694,8 @@ def run_DP3(self, phaseshift: str = None, freqavg: str = None, print(f"DONE: See output --> sub{object.scale}*.ms") else: - print(f"DONE: Output is SUBTRACT_DATA column in input MS") \ No newline at end of file + print(f"DONE: Output is SUBTRACT_DATA column in input MS") + + +if __name__ == "__main__": + main()