diff --git a/auto_selfcal/__main__.py b/auto_selfcal/__main__.py index 42ce51d..cc19510 100644 --- a/auto_selfcal/__main__.py +++ b/auto_selfcal/__main__.py @@ -26,6 +26,7 @@ parser.add_argument('--do_amp_selfcal', default=True) parser.add_argument('--usermask', default={}, type=ast.literal_eval) # require that it is a CRTF region (CASA region format) parser.add_argument('--usermodel', default={}, type=ast.literal_eval) +parser.add_argument('--uniform_solints', default=False) parser.add_argument('--inf_EB_gaincal_combine', default='scan', type=str) # should we get rid of this option? parser.add_argument('--inf_EB_gaintype', default='G', type=str) parser.add_argument('--inf_EB_override', action='store_true') diff --git a/auto_selfcal/auto_selfcal.py b/auto_selfcal/auto_selfcal.py index 4c38643..3eebbcd 100644 --- a/auto_selfcal/auto_selfcal.py +++ b/auto_selfcal/auto_selfcal.py @@ -38,6 +38,7 @@ def auto_selfcal( inf_EB_gaintype='G', inf_EB_override=False, optimize_spw_combine=True, # if False, will not attempt per spw or per baseband solutions for any solint except inf_EB + uniform_solints=False, gaincal_minsnr=2.0, gaincal_unflag_minsnr=5.0, minsnr_to_proceed=2.95, @@ -451,7 +452,7 @@ def auto_selfcal( vis_for_targets[target][band]['vislist'], spectral_average=spectral_average, sort_targets_and_EBs=sort_targets_and_EBs, scale_fov=scale_fov, inf_EB_gaincal_combine=inf_EB_gaincal_combine, inf_EB_gaintype=inf_EB_gaintype, apply_cal_mode_default=apply_cal_mode_default, do_amp_selfcal=do_amp_selfcal, - usermask=usermask, usermodel=usermodel,guess_scan_combine=guess_scan_combine,max_solint=max_solint, + uniform_solints=uniform_solints, usermask=usermask, usermodel=usermodel,guess_scan_combine=guess_scan_combine,max_solint=max_solint, iscalibrator=iscalibrator, do_delay_cal=do_delay_cal, shorter_amp_solints=shorter_amp_solints, imsize=imsize, cell=cell, refant=refant, debug=debug) diff --git a/auto_selfcal/gaincal_wrapper.py b/auto_selfcal/gaincal_wrapper.py index 96401eb..3cf72b6 100644 --- a/auto_selfcal/gaincal_wrapper.py +++ b/auto_selfcal/gaincal_wrapper.py @@ -2,7 +2,7 @@ from .selfcal_helpers import * from .mosaic_helpers import scan_inf_scan_combine -def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, solint_interval, applymode, iteration, +def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, applymode, iteration, gaincal_minsnr, gaincal_unflag_minsnr=5.0, minsnr_to_proceed=3.0, rerank_refants=False, unflag_only_lbants=False, unflag_only_lbants_onlyap=False, calonly_max_flagged=0.0, second_iter_solmode="", unflag_fb_to_prev_solint=False, \ refantmode="flex", mode="selfcal", calibrators="", gaincalibrator_dict={}, allow_gain_interpolation=False,spectral_solution_fraction=0.3, @@ -27,6 +27,8 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so ## Solve gain solutions per MS, target, solint, and band ## + print(selfcal_plan['solmode']) + print(iteration) os.system('rm -rf '+sani_target+'_'+vis+'_'+band+'_'+solint+'_'+str(iteration)+'_'+selfcal_plan['solmode'][iteration]+'*.g') ## Reset the gaincal return dictionaries, in case this is a repeat of the current solution interval. @@ -92,7 +94,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so include_scans = "" print('Include targets: ', include_targets) - if solint == "scan_inf": + if selfcal_plan[vis]['solint_settings'][solint]['sub-name'] == "scan_inf": include_scans = scan_inf_scan_combine(selfcal_library, vis, target, gaincalibrator_dict, guess_scan_combine=guess_scan_combine) else: include_scans = [include_scans] @@ -107,10 +109,10 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so for incl_scan in include_scans: scan_targets = [] for fid in [selfcal_library['sub-fields-fid_map'][vis][fid] for fid in \ - np.intersect1d(selfcal_library['sub-fields-to-gaincal'],list(selfcal_library['sub-fields-fid_map'][vis].keys()))] if incl_scan == '' else \ + np.intersect1d(selfcal_library[vis]['sub-fields-to-gaincal'],list(selfcal_library['sub-fields-fid_map'][vis].keys()))] if incl_scan == '' else \ np.intersect1d(msmd.fieldsforscans(np.array(incl_scan.split(",")).astype(int)), \ [selfcal_library['sub-fields-fid_map'][vis][fid] for fid in \ - numpy.intersect1d(selfcal_library['sub-fields-to-gaincal'],list(selfcal_library['sub-fields-fid_map'][vis].keys()))]): + numpy.intersect1d(selfcal_library[vis]['sub-fields-to-gaincal'],list(selfcal_library['sub-fields-fid_map'][vis].keys()))]): # Note: because of the msmd above getting actual fids from the MS, we just need to append fid below. scan_targets.append(fid) @@ -162,13 +164,14 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so spwselect=','.join(str(spw) for spw in selfcal_library['spws_set'][vis][i].tolist()) else: spwselect=selfcal_library[vis]['spws'] + gaintable_name=sani_target+'_'+vis+'_'+band+'_'+solint+'_'+str(iteration)+'_'+selfcal_plan['solmode'][iteration]+'_'+filename_append+'.g' print('prior to gaincal',gaintable_name, gc_mode) if gc_mode != 'per_bb': gcdict=call_gaincal(vis=vis, caltable=gaintable_name, gaintype=selfcal_plan[vis]['solint_settings'][solint]['gaincal_gaintype'][gc_mode], spw=spwselect, refant=selfcal_library[vis]['refant'], calmode=selfcal_plan['solmode'][iteration], solnorm=solnorm if not do_fallback_calonly else False, - solint=solint_interval.replace('_EB','').replace('_ap','').replace('scan_','').replace('_fb1','').replace('_fb2','').replace('_fb3',''),\ + solint=selfcal_plan[vis]['solint_settings'][solint]['interval'].replace('_EB','').replace('_ap','').replace('scan_','').replace('_fb1','').replace('_fb2','').replace('_fb3',''),\ minsnr=gaincal_minsnr if not do_fallback_calonly else max(gaincal_minsnr,gaincal_unflag_minsnr), minblperant=4,combine=gaincal_combine,\ field=incl_targets,scan=incl_scans,gaintable=gaincal_preapply_gaintable,spwmap=gaincal_spwmap,uvrange=selfcal_library['uvrange'],\ interp=gaincal_interpolate, solmode=gaincal_solmode, refantmode='flex',\ @@ -179,7 +182,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so spwselect_bb=selfcal_library[vis]['baseband'][baseband]['spwstring'] gcdict=call_gaincal(vis=vis, caltable=gaintable_name, gaintype=selfcal_plan[vis]['solint_settings'][solint]['gaincal_gaintype'][gc_mode], spw=spwselect_bb, refant=selfcal_library[vis]['refant'], calmode=selfcal_plan['solmode'][iteration], solnorm=solnorm if not do_fallback_calonly else False, - solint=solint_interval.replace('_EB','').replace('_ap','').replace('scan_','').replace('_fb1','').replace('_fb2','').replace('_fb3',''),\ + solint=selfcal_plan[vis]['solint_settings'][solint]['interval'].replace('_EB','').replace('_ap','').replace('scan_','').replace('_fb1','').replace('_fb2','').replace('_fb3',''),\ minsnr=gaincal_minsnr if not do_fallback_calonly else max(gaincal_minsnr,gaincal_unflag_minsnr), minblperant=4,combine=gaincal_combine,\ field=incl_targets,scan=incl_scans,gaintable=gaincal_preapply_gaintable,spwmap=gaincal_spwmap,uvrange=selfcal_library['uvrange'],\ interp=gaincal_interpolate, solmode=gaincal_solmode, refantmode='flex',\ @@ -195,7 +198,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so gaintable_prefix=sani_target+'_'+vis+'_'+band+'_' # assume that if there is only one mode to attempt, that it is combinespw and don't bother checking. - if 'delay' not in solint: + if 'd' not in solint: if len(selfcal_plan[vis]['solint_settings'][solint]['modes_to_attempt']) > 1: get_gaincalmode_flagging_stats(selfcal_library,selfcal_plan,vis,gaintable_prefix,solint) preferred_mode,fallback,spwmap,spwmapping_for_applycal = \ @@ -244,7 +247,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so remove_modes(selfcal_plan,vis,current_solint_index) - for fid in np.intersect1d(selfcal_library['sub-fields-to-selfcal'],list(selfcal_library['sub-fields-fid_map'][vis].keys())): + for fid in np.intersect1d(selfcal_library[vis]['sub-fields-to-selfcal'],list(selfcal_library['sub-fields-fid_map'][vis].keys())): selfcal_library[fid][vis][solint]['final_mode']=preferred_mode+'' selfcal_library[fid][vis][solint]['spwmap']=applycal_spwmap.copy() selfcal_library[fid][vis][solint]['gaincal_combine']=selfcal_plan[vis]['solint_settings'][solint]['gaincal_combine'][preferred_mode]+'' @@ -263,7 +266,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so else: # this else is for ap selfcal os.system('rm -rf temp*.g') - for fid in np.intersect1d(selfcal_library['sub-fields-to-selfcal'],list(selfcal_library['sub-fields-fid_map'][vis].keys())): + for fid in np.intersect1d(selfcal_library[vis]['sub-fields-to-selfcal'],list(selfcal_library['sub-fields-fid_map'][vis].keys())): gaincal_spwmap=[] gaincal_preapply_gaintable=[] gaincal_interpolate=[] @@ -271,6 +274,8 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so applycal_interpolate=[] applycal_spwmap=[] for j in range(current_solint_index): + if not selfcal_plan['solints'][j] in selfcal_plan[vis]['solint_settings']: + continue if selfcal_plan[vis]['solint_settings'][selfcal_plan['solints'][j]]['preapply_this_gaintable']: #allow amplitude and phase to be pre-applied # and selfcal_plan[vis]['solint_settings'][selfcal_plan['solints'][j]]['solmode']=='p': gaincal_preapply_gaintable.append(selfcal_plan[vis]['solint_settings'][selfcal_plan['solints'][j]]['accepted_gaintable']) @@ -306,11 +311,12 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so spwselect=','.join(str(spw) for spw in selfcal_library['spws_set'][vis][i].tolist()) else: spwselect=selfcal_library[vis]['spws'] + gaintable_name='temp_'+solint+'_'+str(iteration)+'_'+selfcal_plan['solmode'][iteration]+'_'+filename_append+'.g' if gc_mode != 'per_bb': gcdict=call_gaincal(vis=vis, caltable=gaintable_name, gaintype=selfcal_plan[vis]['solint_settings'][solint]['gaincal_gaintype'][gc_mode], spw=spwselect, refant=selfcal_library[vis]['refant'], calmode=selfcal_plan['solmode'][iteration], solnorm=solnorm if not do_fallback_calonly else False, - solint=solint_interval.replace('_EB','').replace('_ap','').replace('scan_',''),\ + solint=selfcal_plan[vis]['solint_settings'][solint]['interval'].replace('_EB','').replace('_ap','').replace('scan_',''),\ minsnr=gaincal_minsnr if not do_fallback_calonly else max(gaincal_minsnr,gaincal_unflag_minsnr), minblperant=4,combine=gaincal_combine,\ field=str(selfcal_library['sub-fields-fid_map'][vis][fid]),gaintable=gaincal_preapply_gaintable,spwmap=gaincal_spwmap,uvrange=selfcal_library['uvrange'],\ interp=gaincal_interpolate, solmode=gaincal_solmode, refantmode='flex',\ @@ -321,7 +327,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so spwselect_bb=selfcal_library[vis]['baseband'][baseband]['spwstring'] gcdict=call_gaincal(vis=vis, caltable=gaintable_name, gaintype=selfcal_plan[vis]['solint_settings'][solint]['gaincal_gaintype'][gc_mode], spw=spwselect_bb, refant=selfcal_library[vis]['refant'], calmode=selfcal_plan['solmode'][iteration], solnorm=solnorm if not do_fallback_calonly else False, - solint=solint_interval.replace('_EB','').replace('_ap','').replace('scan_',''),\ + solint=selfcal_plan[vis]['solint_settings'][solint]['interval'].replace('_EB','').replace('_ap','').replace('scan_',''),\ minsnr=gaincal_minsnr if not do_fallback_calonly else max(gaincal_minsnr,gaincal_unflag_minsnr), minblperant=4,combine=gaincal_combine,\ field=str(selfcal_library['sub-fields-fid_map'][vis][fid]),gaintable=gaincal_preapply_gaintable,spwmap=gaincal_spwmap,uvrange=selfcal_library['uvrange'],\ interp=gaincal_interpolate, solmode=gaincal_solmode, refantmode='flex',\ @@ -337,7 +343,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so select_best_gaincal_mode(selfcal_library,selfcal_plan,vis,gaintable_prefix,solint,spectral_solution_fraction,minsnr_to_proceed) else: preferred_mode=selfcal_plan[vis]['solint_settings'][solint]['modes_to_attempt'][0] - if 'delay' not in solint: + if 'd' not in solint: get_gaincalmode_flagging_stats(selfcal_library,selfcal_plan,vis,gaintable_prefix,solint) fallback='' spwmapping_for_applycal=[] @@ -367,7 +373,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so remove_modes(selfcal_plan,vis,current_solint_index) - for fid in np.intersect1d(selfcal_library['sub-fields-to-selfcal'],list(selfcal_library['sub-fields-fid_map'][vis].keys())): + for fid in np.intersect1d(selfcal_library[vis]['sub-fields-to-selfcal'],list(selfcal_library['sub-fields-fid_map'][vis].keys())): selfcal_library[fid][vis][solint]['final_mode']=preferred_mode+'' selfcal_library[fid][vis][solint]['spwmap']=applycal_spwmap.copy() selfcal_library[fid][vis][solint]['gaincal_combine']=selfcal_plan[vis]['solint_settings'][solint]['gaincal_combine'][preferred_mode]+'' @@ -442,7 +448,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so refant=selfcal_library[vis]["refant"], refantmode=refantmode if 'inf_EB' not in solint else 'flex') selfcal_library[vis][solint]['fallback']=fallback+'' - for fid in np.intersect1d(selfcal_library['sub-fields-to-selfcal'],list(selfcal_library['sub-fields-fid_map'][vis].keys())): + for fid in np.intersect1d(selfcal_library[vis]['sub-fields-to-selfcal'],list(selfcal_library['sub-fields-fid_map'][vis].keys())): selfcal_library[fid][vis][solint]['fallback']=fallback+'' # If iteration two, try restricting to just the antennas with enough unflagged data. @@ -486,7 +492,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so if (solint != "inf_EB" and not allow_gain_interpolation) or (allow_gain_interpolation and "inf" not in solint): # If a given field has > 25% of its solutions flagged then just flag the whole field because it will have too much # interpolation. - if solint == "scan_inf": + if selfcal_plan[vis]['solint_settings'][solint]['sub-name'] == "scan_inf": max_n_solutions = max([(scans == scan).sum() for scan in np.unique(scans)]) for scan in np.unique(scans): scan_n_solutions = (flags[0,0,scans == scan] == False).sum() @@ -509,7 +515,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so def generate_settings_for_combinespw_fallback(selfcal_library, selfcal_plan, target, band, vis, solint, iteration): sani_target=sanitize_string(target) current_solint_index=selfcal_plan['solints'].index(solint) - if selfcal_library['telescope'] == 'VLBA' or 'delay' in solint: # use per_bb in place of combinespw for VLBA fall back' + if selfcal_library['telescope'] == 'VLBA' or 'd' in solint: # use per_bb in place of combinespw for VLBA fall back' preferred_mode='per_bb' else: preferred_mode='combinespw' @@ -544,7 +550,7 @@ def generate_settings_for_combinespw_fallback(selfcal_library, selfcal_plan, tar selfcal_library[vis][solint]['applycal_interpolate']=applycal_interpolate selfcal_library[vis][solint]['solmode']=selfcal_plan['solmode'][iteration]+'' selfcal_library[vis][solint]['gaincal_combine']=selfcal_plan[vis]['solint_settings'][solint]['gaincal_combine'][preferred_mode]+'' - for fid in np.intersect1d(selfcal_library['sub-fields-to-selfcal'],list(selfcal_library['sub-fields-fid_map'][vis].keys())): + for fid in np.intersect1d(selfcal_library[vis]['sub-fields-to-selfcal'],list(selfcal_library['sub-fields-fid_map'][vis].keys())): selfcal_library[fid][vis][solint]['final_mode']=preferred_mode+'_fallback' selfcal_library[fid][vis][solint]['spwmap']=applycal_spwmap selfcal_library[fid][vis][solint]['gaincal_combine']=selfcal_plan[vis]['solint_settings'][solint]['gaincal_combine'][preferred_mode]+'' diff --git a/auto_selfcal/image_analysis_helpers.py b/auto_selfcal/image_analysis_helpers.py index 20b0a77..a469093 100644 --- a/auto_selfcal/image_analysis_helpers.py +++ b/auto_selfcal/image_analysis_helpers.py @@ -12,7 +12,12 @@ def get_image_stats(image, mask, backup_mask, selfcal_library, use_nfmask, solin else: SNR_NF, RMS_NF = SNR, RMS - for vis in selfcal_library['vislist']: + if suffix in ['dirty','orig','initial','final']: + vislist = selfcal_library['vislist'] + else: + vislist = selfcal_library['vislist-to-gaincal'] + + for vis in vislist: if suffix in ['dirty','orig','initial','final']: if spw == 'all': update_dict = selfcal_library diff --git a/auto_selfcal/mosaic_helpers.py b/auto_selfcal/mosaic_helpers.py index 1db762f..06715f2 100644 --- a/auto_selfcal/mosaic_helpers.py +++ b/auto_selfcal/mosaic_helpers.py @@ -1,14 +1,14 @@ import numpy as np from .selfcal_helpers import * -def evaluate_subfields_to_gaincal(selfcal_library, target, band, solint, iteration, solmode, solints, selfcal_plan, +def evaluate_subfields_to_gaincal(vis, selfcal_library, target, band, solint, iteration, solmode, solints, selfcal_plan, minsnr_to_proceed, allow_gain_interpolation=False,): sani_target=sanitize_string(target) # Fields that don't have any mask in the primary beam should be removed from consideration, as their models are likely bad. new_fields_to_selfcal = [] - for fid in selfcal_library['sub-fields-to-selfcal']: + for fid in selfcal_library[vis]['sub-fields-to-selfcal']: os.system('rm -rf test*.mask') tmp_SNR_NF,tmp_RMS_NF=estimate_near_field_SNR(sani_target+'_field_'+str(fid)+'_'+band+'_'+solint+'_'+str(iteration)+'.image.tt0', \ las=selfcal_library['LAS'], mosaic_sub_field=True, save_near_field_mask=False) @@ -41,8 +41,8 @@ def evaluate_subfields_to_gaincal(selfcal_library, target, band, solint, iterati if not checkmask(sani_target+'_field_'+str(fid)+'_'+band+'_'+solint+'_'+str(iteration)+'.image.tt0'): print("Removing field "+str(fid)+" from gaincal because there is no signal within the primary beam.") skip_reason = "No signal" - elif selfcal_plan[fid]['solint_snr_per_field'][solints[iteration]] < minsnr_to_proceed and \ - solint not in ['inf_EB','scan_inf']: + elif selfcal_plan[fid][vis]['solint_snr_per_field'][solints[iteration]] < minsnr_to_proceed and \ + selfcal_plan[vis]['solint_settings'][solint]['sub-name'] not in ['inf_EB','scan_inf']: print("Removing field "+str(fid)+" from gaincal because the estimated solint snr is too low.") skip_reason = "Estimated SNR" elif updated_intflux > selfcal_library['flux_threshold'] * original_intflux: @@ -53,21 +53,20 @@ def evaluate_subfields_to_gaincal(selfcal_library, target, band, solint, iterati new_fields_to_selfcal.append(fid) if fid not in new_fields_to_selfcal and solint != "inf_EB" and not allow_gain_interpolation: - for vis in selfcal_library[fid]['vislist']: - #selfcal_library[fid][vis][solint]['interpolated_gains'] = True - #selfcal_library[fid]['Stop_Reason'] = "Gaincal solutions would be interpolated" - selfcal_library[fid]['Stop_Reason'] = skip_reason - selfcal_library[fid][vis][solint]['Pass'] = "None" - selfcal_library[fid][vis][solint]['Fail_Reason'] = skip_reason + #selfcal_library[fid][vis][solint]['interpolated_gains'] = True + #selfcal_library[fid]['Stop_Reason'] = "Gaincal solutions would be interpolated" + selfcal_library[fid][vis]['Stop_Reason'] = skip_reason + selfcal_library[fid][vis][solint]['Pass'] = "None" + selfcal_library[fid][vis][solint]['Fail_Reason'] = skip_reason return new_fields_to_selfcal -def evaluate_subfields_after_gaincal(selfcal_library, target, band, solint, iteration, solmode, allow_gain_interpolation=False): +def evaluate_subfields_after_gaincal(vis, selfcal_library, selfcal_plan, target, band, solint, iteration, solmode, allow_gain_interpolation=False): - new_fields_to_selfcal = selfcal_library['sub-fields-to-selfcal'].copy() + new_fields_to_selfcal = selfcal_library[vis]['sub-fields-to-selfcal'].copy() sani_target=sanitize_string(target) @@ -75,23 +74,23 @@ def evaluate_subfields_after_gaincal(selfcal_library, target, band, solint, iter (allow_gain_interpolation and "inf" not in solint)): # With gaincal done and bad fields removed from gain tables if necessary, check whether any fields should no longer be selfcal'd # because they have too much interpolation. - for vis in selfcal_library['vislist']: + #for vis in selfcal_library['vislist']: + if True: ## If an EB had no fields to gaincal on, remove all fields in that EB from being selfcal'd as there is no calibration available ## in this EB. - if np.intersect1d(selfcal_library['sub-fields-to-gaincal'],\ + if np.intersect1d(selfcal_library[vis]['sub-fields-to-gaincal'],\ list(selfcal_library['sub-fields-fid_map'][vis].keys())).size == 0: for fid in np.intersect1d(new_fields_to_selfcal,list(selfcal_library['sub-fields-fid_map'][vis].keys())): new_fields_to_selfcal.remove(fid) - selfcal_library[fid]['Stop_Reason'] = 'No viable calibrator fields in at least 1 EB' - for v in selfcal_library[fid]['vislist']: - selfcal_library[fid][v][solint]['Pass'] = 'None' - if 'Fail_Reason' in selfcal_library[fid][v][solint]: - selfcal_library[fid][v][solint]['Fail_Reason'] += '; ' - else: - selfcal_library[fid][v][solint]['Fail_Reason'] = '' - selfcal_library[fid][v][solint]['Fail_Reason'] += 'No viable fields' - continue + selfcal_library[fid][vis]['Stop_Reason'] = 'No viable calibrator fields in at least 1 EB' + selfcal_library[fid][vis][solint]['Pass'] = 'None' + if 'Fail_Reason' in selfcal_library[fid][vis][solint]: + selfcal_library[fid][vis][solint]['Fail_Reason'] += '; ' + else: + selfcal_library[fid][vis][solint]['Fail_Reason'] = '' + selfcal_library[fid][vis][solint]['Fail_Reason'] += 'No viable fields' + return new_fields_to_selfcal ## NEXT TO DO: check % of flagged solutions - DONE, see above ## After that enable option for interpolation through inf - DONE tb.open(sani_target+'_'+vis+'_'+band+'_'+solint+'_'+str(iteration)+'_'+solmode[iteration]+'_'+ @@ -100,7 +99,7 @@ def evaluate_subfields_after_gaincal(selfcal_library, target, band, solint, iter scans = tb.getcol("SCAN_NUMBER") for fid in np.intersect1d(new_fields_to_selfcal,list(selfcal_library['sub-fields-fid_map'][vis].keys())): - if solint == "scan_inf": + if selfcal_plan[vis]['solint_settings'][solint]['sub-name'] == "scan_inf": msmd.open(vis) cals_for_scan = [] total_cals_for_scan = [] @@ -125,29 +124,29 @@ def evaluate_subfields_after_gaincal(selfcal_library, target, band, solint, iter new_fields_to_selfcal.remove(fid) if fid not in new_fields_to_selfcal: - # We need to update all the EBs, not just the one that failed. - for v in selfcal_library[fid]['vislist']: - selfcal_library[fid][v][solint]['Pass'] = 'None' - if allow_gain_interpolation: - selfcal_library[fid][v][solint]['Fail_Reason'] = 'Interpolation beyond inf' - else: - selfcal_library[fid][v][solint]['Fail_Reason'] = 'Bad gaincal solutions' + selfcal_library[fid][vis][solint]['Pass'] = 'None' + if allow_gain_interpolation: + selfcal_library[fid][vis][solint]['Fail_Reason'] = 'Interpolation beyond inf' + else: + selfcal_library[fid][vis][solint]['Fail_Reason'] = 'Bad gaincal solutions' tb.close() elif selfcal_library['obstype'] == 'mosaic' and solint == "inf_EB": ## If an EB had no fields to gaincal on, remove all fields in that EB from being selfcal'd as there is no calibration available ## in this EB. - for vis in selfcal_library['vislist']: - if np.intersect1d(selfcal_library['sub-fields-to-gaincal'],\ + #for vis in selfcal_library['vislist']: + if True: + if np.intersect1d(selfcal_library[vis]['sub-fields-to-gaincal'],\ list(selfcal_library['sub-fields-fid_map'][vis].keys())).size == 0: for fid in np.intersect1d(new_fields_to_selfcal,list(selfcal_library['sub-fields-fid_map'][vis].keys())): new_fields_to_selfcal.remove(fid) - selfcal_library[fid]['Stop_Reason'] = 'No viable calibrator fields for inf_EB in at least 1 EB' - for v in selfcal_library[fid]['vislist']: - selfcal_library[fid][v][solint]['Pass'] = 'None' - selfcal_library[fid][v][solint]['Fail_Reason'] = 'No viable inf_EB fields' + selfcal_library[fid][vis]['Stop_Reason'] = 'No viable calibrator fields for inf_EB in at least 1 EB' + selfcal_library[fid][vis][solint]['Pass'] = 'None' + selfcal_library[fid][vis][solint]['Fail_Reason'] = 'No viable inf_EB fields' + + print("new_fields_to_selfcal", new_fields_to_selfcal) return new_fields_to_selfcal diff --git a/auto_selfcal/original_ms_helpers.py b/auto_selfcal/original_ms_helpers.py index f41bb32..7a56378 100644 --- a/auto_selfcal/original_ms_helpers.py +++ b/auto_selfcal/original_ms_helpers.py @@ -75,7 +75,7 @@ def applycal_to_orig_MSes(selfcal_library='selfcal_library.pickle', write_only=T for band in selfcal_library[target].keys(): if selfcal_library[target][band]['SC_success']: for vis in selfcal_library[target][band]['vislist']: - solint=selfcal_library[target][band]['final_solint'] + solint=selfcal_library[target][band][vis]['final_solint'] iteration=selfcal_library[target][band][vis][solint]['iteration'] # Write the line to the command log diff --git a/auto_selfcal/prepare_cocal.py b/auto_selfcal/prepare_cocal.py index 0d61498..8871750 100644 --- a/auto_selfcal/prepare_cocal.py +++ b/auto_selfcal/prepare_cocal.py @@ -55,7 +55,12 @@ def prepare_cocal(selfcal_library, selfcal_plan, inf_EB_gaincal_combine, inf_EB_ # Update the relevant lists if we are going to do a fallback mode. selfcal_plan[target][band]['solints'] += ["inf_EB_fb","inf_fb1","inf_fb2","inf_fb3"] selfcal_plan[target][band]['solmode'] += ["p","p","p","p"] - selfcal_plan[target][band]['solint_interval'] += ["inf","inf","inf","inf"] + + for vis in selfcal_library[target][band]['vislist']: + for solint_name in ["inf_EB_fb","inf_fb1","inf_fb2","inf_fb3"]: + selfcal_plan[target][band][vis]['solint_settings'][solint_name] = {} + selfcal_plan[target][band][vis]['solint_settings'][solint_name]['interval'] = solint_name + selfcal_plan[target][band][vis]['solint_settings'][solint_name]['sub-name'] = solint_name plan_selfcal_per_solint(selfcal_library, selfcal_plan, @@ -83,7 +88,6 @@ def prepare_cocal(selfcal_library, selfcal_plan, inf_EB_gaincal_combine, inf_EB_ for cal_target in inf_fields[band]: selfcal_plan[target][band][vis]['solint_settings']["inf_fb3"]["preapply_gaintable_dict"][cal_target] = selfcal_plan[cal_target][band][vis.replace(sanitize_string(target), sanitize_string(cal_target))]['solint_settings']['inf_EB']['accepted_gaintable'] - selfcal_plan[target][band]['gaincal_combine'] += [selfcal_plan[target][band]['gaincal_combine'][0], selfcal_plan[target][band]['gaincal_combine'][1], selfcal_plan[target][band]['gaincal_combine'][1], selfcal_plan[target][band]['gaincal_combine'][1]] selfcal_plan[target][band]['applycal_mode'] += [selfcal_plan[target][band]['applycal_mode'][0], selfcal_plan[target][band]['applycal_mode'][1], selfcal_plan[target][band]['applycal_mode'][1], selfcal_plan[target][band]['applycal_mode'][1]] calibrators[band] = [inf_EB_fields[band], inf_fields[band], inf_fields[band], inf_fields[band]] diff --git a/auto_selfcal/prepare_selfcal.py b/auto_selfcal/prepare_selfcal.py index bb24888..54e3e5f 100644 --- a/auto_selfcal/prepare_selfcal.py +++ b/auto_selfcal/prepare_selfcal.py @@ -11,6 +11,7 @@ def prepare_selfcal(all_targets, bands, bands_for_targets, vislist, inf_EB_gaintype='G', apply_cal_mode_default='calflag', do_amp_selfcal=True, + uniform_solints=False, usermask={}, usermodel={}, max_solint=4500.0, @@ -151,6 +152,8 @@ def prepare_selfcal(all_targets, bands, bands_for_targets, vislist, selfcal_library[target][band]['sub-fields'] = list(range(len(all_phasecenters))) selfcal_library[target][band]['sub-fields-to-selfcal'] = list(range(len(all_phasecenters))) + for vis in vislist: + selfcal_library[target][band][vis]['sub-fields-to-selfcal'] = list(range(len(all_phasecenters))) selfcal_library[target][band]['sub-fields-phasecenters'] = dict(zip(selfcal_library[target][band]['sub-fields'], all_phasecenters)) # Now we can start to create a sub-field selfcal_library entry for each sub-field. @@ -160,6 +163,8 @@ def prepare_selfcal(all_targets, bands, bands_for_targets, vislist, for vis in vislist: if not fid in selfcal_library[target][band]['sub-fields-fid_map'][vis]: + # If a sub-field is not in an EB, it shouldn't be considered for selfcal for that EB + selfcal_library[target][band][vis]['sub-fields-to-selfcal'].remove(fid) continue selfcal_library[target][band][fid][vis] = {} @@ -502,16 +507,78 @@ def default(self, obj): selfcal_plan[target] = {} for band in selfcal_library[target]: - selfcal_plan[target][band] = {} if band in selfcal_library[target]: - selfcal_plan[target][band]['solints'],selfcal_plan[target][band]['integration_time'],selfcal_plan[target][band]['gaincal_combine'], \ - selfcal_plan[target][band]['solmode'],selfcal_plan[target][band]['solint_interval']=get_solints_simple(selfcal_library[target][band]['vislist'],scantimesdict[band],\ - scannfieldsdict[band],scanstartsdict[band],scanendsdict[band],integrationtimesdict[band],\ - inf_EB_gaincal_combine,do_amp_selfcal=do_amp_selfcal,mosaic=selfcal_library[target][band]['obstype'] == 'mosaic',do_scan_inf=do_scan_inf,\ - max_solint=max_solint,iscalibrator=iscalibrator,shorter_amp_solints=shorter_amp_solints,do_delay_cal=do_delay_cal,n_solints=n_solints) + selfcal_plan[target][band] = {} + selfcal_plan[target][band]['solints'] = [] + selfcal_plan[target][band]['solmode'] = [] + + if uniform_solints: + solints,tmp_integration_time,tmp_gaincal_combine, \ + tmp_solmodes,tmp_solint_interval=get_solints_simple(vislist,scantimesdict[band], + scannfieldsdict[band],scanstartsdict[band],scanendsdict[band],integrationtimesdict[band],\ + inf_EB_gaincal_combine,do_amp_selfcal=do_amp_selfcal,mosaic=selfcal_library[target][band]['obstype'] == 'mosaic',do_scan_inf=do_scan_inf,\ + max_solint=max_solint,iscalibrator=iscalibrator,shorter_amp_solints=shorter_amp_solints,do_delay_cal=do_delay_cal,n_solints=n_solints) + + for vis in selfcal_library[target][band]['vislist']: + selfcal_plan[target][band][vis] = {} + + if not uniform_solints: + solints,selfcal_plan[target][band][vis]['integration_time'],selfcal_plan[target][band][vis]['gaincal_combine'], \ + tmp_solmodes,tmp_solint_interval=get_solints_simple([vis],scantimesdict[band], + scannfieldsdict[band],scanstartsdict[band],scanendsdict[band],integrationtimesdict[band],\ + inf_EB_gaincal_combine,do_amp_selfcal=do_amp_selfcal,mosaic=selfcal_library[target][band]['obstype'] == 'mosaic',do_scan_inf=do_scan_inf,\ + max_solint=max_solint,iscalibrator=iscalibrator,shorter_amp_solints=shorter_amp_solints,do_delay_cal=do_delay_cal,n_solints=n_solints) + else: + selfcal_plan[target][band][vis]['integration_time'] = tmp_integration_time + selfcal_plan[target][band][vis]['gaincal_combine'] = tmp_gaincal_combine + + selfcal_plan[target][band][vis]['solint_settings']={} + + subscan_count = 0 + for solint, solint_interval in zip(solints, tmp_solint_interval): + """ + if 'inf' in solint or 'int' in solint or 'ap' in solint: + solint_name = solint + else: + solint_name = 'subscan'+str(subscan_count) + subscan_count += 1 + """ + if solint == "inf_EB" or solint == "inf_EB_delay": + solint_name = solint + else: + if 'ap' in solint: + solint_name = 'ap' + elif 'delay' in solint: + solint_name = 'd' + else: + solint_name = 'p' + solint_name += str(subscan_count) + if solint == 'int': + subscan_count = 0 + else: + subscan_count += 1 + + if solint_name not in selfcal_plan[target][band]['solints']: + if 'ap' not in solint_name and 'ap' in selfcal_plan[target][band]['solmode']: + insert_loc = np.where(np.array(selfcal_plan[target][band]['solmode']) == 'ap')[0][0] + else: + insert_loc = len(selfcal_plan[target][band]['solmode']) + + selfcal_plan[target][band]['solints'].insert(insert_loc,solint_name) + if 'ap' in solint_name: + selfcal_plan[target][band]['solmode'].insert(insert_loc,'ap') + else: + selfcal_plan[target][band]['solmode'].insert(insert_loc,'p') + + selfcal_plan[target][band][vis]['solint_settings'][solint_name]={} + selfcal_plan[target][band][vis]['solint_settings'][solint_name]['interval'] = solint_interval + selfcal_plan[target][band][vis]['solint_settings'][solint_name]['sub-name'] = solint + + selfcal_plan[target][band]['applycal_mode']=[apply_cal_mode_default]*len(selfcal_plan[target][band]['solints']) + print(band,target,selfcal_plan[target][band]['solints']) - print(band,target,selfcal_plan[target][band]['solint_interval']) - selfcal_plan[target][band]['applycal_mode']=[apply_cal_mode_default]*len(selfcal_plan[target][band]['solints']) + for vis in vislist: + print(vis,[selfcal_plan[target][band][vis]['solint_settings'][solint]['interval'] for solint in selfcal_plan[target][band]['solints'] if solint in selfcal_plan[target][band][vis]['solint_settings']]) ## ## estimate per scan/EB S/N using time on source and median scan times @@ -520,7 +587,6 @@ def default(self, obj): for target in selfcal_plan: for band in selfcal_plan[target]: for vis in selfcal_library[target][band]['vislist']: - selfcal_plan[target][band][vis] = {} selfcal_plan[target][band][vis]['inf_EB_gaincal_combine']=inf_EB_gaincal_combine #'scan' if selfcal_library[target][band]['obstype']=='mosaic': selfcal_plan[target][band][vis]['inf_EB_gaincal_combine']+=',field' @@ -562,13 +628,15 @@ def plan_selfcal_per_solint(selfcal_library, selfcal_plan, optimize_spw_combine= if solints is not None: use_solints = solints else: - selfcal_plan[target][band][vis]['solint_settings']={} + #selfcal_plan[target][band][vis]['solint_settings']={} use_solints = selfcal_plan[target][band]['solints'] for solint in use_solints: + if solint not in selfcal_plan[target][band][vis]['solint_settings']: + continue gaincal_combine='' filename_append='' - selfcal_plan[target][band][vis]['solint_settings'][solint]={} + #selfcal_plan[target][band][vis]['solint_settings'][solint]={} selfcal_plan[target][band][vis]['solint_settings'][solint]['preapply_this_gaintable']=False selfcal_plan[target][band][vis]['solint_settings'][solint]['gaincal_preapply_gaintable']=[] selfcal_plan[target][band][vis]['solint_settings'][solint]['gaincal_spwmap']=[] @@ -580,7 +648,7 @@ def plan_selfcal_per_solint(selfcal_library, selfcal_plan, optimize_spw_combine= selfcal_plan[target][band][vis]['solint_settings'][solint]['applycal_gaintable']=[] selfcal_plan[target][band][vis]['solint_settings'][solint]['applycal_spwmap']=[] selfcal_plan[target][band][vis]['solint_settings'][solint]['spwmap_for_mode']={} - if 'delay' not in solint: + if 'd' not in solint: selfcal_plan[target][band][vis]['solint_settings'][solint]['applycal_interpolate']=applycal_interp else: selfcal_plan[target][band][vis]['solint_settings'][solint]['applycal_interpolate']='linear' @@ -589,11 +657,11 @@ def plan_selfcal_per_solint(selfcal_library, selfcal_plan, optimize_spw_combine= selfcal_plan[target][band][vis]['solint_settings'][solint]['modes_to_attempt']=[] print('Nspws: {}, spws per BB: {}, basebands: {}'.format(nspws,maxspws_per_bb,n_basebands)) - if selfcal_library[target][band]['telescope'] == 'VLBA' and 'delay' in solint and maxspws_per_bb > 1.0: + if selfcal_library[target][band]['telescope'] == 'VLBA' and 'd' in solint and maxspws_per_bb > 1.0: selfcal_plan[target][band][vis]['solint_settings'][solint]['modes_to_attempt'].append('per_bb') - if 'delay' not in solint and nspws > 1.0: + if 'd' not in solint and nspws > 1.0: selfcal_plan[target][band][vis]['solint_settings'][solint]['modes_to_attempt'].append('combinespw') - if 'delay' in solint and n_basebands > 1: + if 'd' in solint and n_basebands > 1: selfcal_plan[target][band][vis]['solint_settings'][solint]['modes_to_attempt'].append('per_bb') if "inf_EB" in solint: selfcal_plan[target][band][vis]['solint_settings'][solint]['modes_to_attempt'].append('combinespwpol') @@ -601,8 +669,8 @@ def plan_selfcal_per_solint(selfcal_library, selfcal_plan, optimize_spw_combine= selfcal_plan[target][band][vis]['solint_settings'][solint]['preapply_this_gaintable']=True print('solint',solint,'N basebands: ',n_basebands, 'modes to attempt: ',selfcal_plan[target][band][vis]['solint_settings'][solint]['modes_to_attempt']) if 'spw' not in selfcal_plan[target][band][vis]['inf_EB_gaincal_combine'] and 'fb' not in solint: - min_SNR_spw=get_min_SNR_spw(selfcal_plan[target][band]['solint_snr_per_spw'][solint]) - min_SNR_bb=get_min_SNR_spw(selfcal_plan[target][band]['solint_snr_per_bb'][solint]) + min_SNR_spw=get_min_SNR_spw(selfcal_plan[target][band][vis]['solint_snr_per_spw'][solint]) + min_SNR_bb=get_min_SNR_spw(selfcal_plan[target][band][vis]['solint_snr_per_bb'][solint]) if min_SNR_spw > 2.0: selfcal_plan[target][band][vis]['solint_settings'][solint]['modes_to_attempt'].append('per_spw') #selfcal_plan[target][band][vis]['solint_settings'][solint]['preapply_this_gaintable']=True # leave default to off and have it decide after eval @@ -614,7 +682,7 @@ def plan_selfcal_per_solint(selfcal_library, selfcal_plan, optimize_spw_combine= selfcal_plan[target][band][vis]['solint_settings'][solint]['modes_to_attempt'].append('per_bb') #selfcal_plan[target][band][vis]['solint_settings'][solint]['preapply_this_gaintable']=True # leave default to off and have it decide after eval - if '_ap' in solint: + if 'ap' in solint: selfcal_plan[target][band][vis]['solint_settings'][solint]['solmode']='ap' else: selfcal_plan[target][band][vis]['solint_settings'][solint]['solmode']='p' @@ -640,11 +708,11 @@ def plan_selfcal_per_solint(selfcal_library, selfcal_plan, optimize_spw_combine= gaincal_combine='spw' filename_append='per_bb' selfcal_plan[target][band][vis]['solint_settings'][solint]['spwmap_for_mode']['per_bb']=selfcal_library[target][band][vis]['baseband_spwmap'] - if solint in ['inf_EB','inf_EB_delay','scan_inf','300s_ap','inf_EB_fb']: + if selfcal_plan[target][band][vis]['solint_settings'][solint]['sub-name'] in ['inf_EB','inf_EB_delay','scan_inf','300s_ap','inf_EB_fb']: if gaincal_combine!='': gaincal_combine+=',' gaincal_combine+='scan' - if (solint in ['inf_EB','inf_EB_delay','scan_inf'] and selfcal_library[target][band]['obstype'] == 'mosaic') or \ + if (selfcal_plan[target][band][vis]['solint_settings'][solint]['sub-name'] in ['inf_EB','inf_EB_delay','scan_inf'] and selfcal_library[target][band]['obstype'] == 'mosaic') or \ solint == "inf_EB_fb": gaincal_combine+=',field' selfcal_plan[target][band][vis]['solint_settings'][solint]['gaincal_combine'][mode]=gaincal_combine @@ -659,7 +727,7 @@ def plan_selfcal_per_solint(selfcal_library, selfcal_plan, optimize_spw_combine= selfcal_plan[target][band][vis]['solint_settings'][solint]['gaincal_gaintype'][mode]='T' if selfcal_library[target][band]['telescope'] == 'VLBA': # do G by default for VLBA selfcal_plan[target][band][vis]['solint_settings'][solint]['gaincal_gaintype'][mode]='G' - if '_delay' in solint : # delay solints should always use K and be preapplied + if 'd' in solint : # delay solints should always use K and be preapplied selfcal_plan[target][band][vis]['solint_settings'][solint]['gaincal_gaintype'][mode]='K' #consider if improvement from delay is neglegible or harms to continume selfcal but then turning off the preapply selfcal_plan[target][band][vis]['solint_settings'][solint]['preapply_this_gaintable']=True diff --git a/auto_selfcal/run_selfcal.py b/auto_selfcal/run_selfcal.py index 93e962c..97a6b84 100644 --- a/auto_selfcal/run_selfcal.py +++ b/auto_selfcal/run_selfcal.py @@ -44,11 +44,9 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ else: iteration = 0 - vislist=selfcal_library['vislist'].copy() - if mode == "cocal": # Check whether there are suitable calibrators, otherwise skip this target/band. - include_targets, include_scans = triage_calibrators(vislist[0], target, band, calibrators[band][0]) + include_targets, include_scans = triage_calibrators(selfcal_library['vislist'][0], target, band, calibrators[band][0]) print('Co-calibrators: ',include_targets) print('Co-calibrator scans: ',include_scans) if include_targets == "": @@ -56,7 +54,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ selfcal_library['Stop_Reason'] += '; No suitable co-calibrators' return else: - for vis in vislist: + for vis in selfcal_library['vislist']: clearcal(vis=vis,addmodel=True) os.system('mv '+vis+' '+vis.replace('.ms','_orig.ms')) os.system('mv '+vis+'.flagversions '+vis.replace('.ms','_orig.ms.flagversions')) @@ -80,11 +78,17 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ do_fallback_combinespw=False do_fallback_calonly=False print('Starting selfcal procedure on: '+target+' '+band) + print('selfcal_library["sub-fields-to-selfcal"] = ', selfcal_library['sub-fields-to-selfcal']) + for vis in selfcal_library['vislist']: + print('selfcal_library[vis]["sub-fields-to-selfcal"] = ', selfcal_library[vis]['sub-fields-to-selfcal']) while iteration < len(selfcal_plan['solints']): - print(selfcal_plan['solints']) - print(selfcal_plan['solint_interval'][iteration]) - print("Solving for solint="+selfcal_plan['solints'][iteration]+' with interval '+selfcal_plan['solint_interval'][iteration]) - + vislist=[vis for vis in selfcal_library['vislist'] if selfcal_plan['solints'][iteration] in selfcal_plan[vis]['solint_settings']] + + print("Solving for solint="+selfcal_plan['solints'][iteration]+' with intervals:') + for vis in vislist: + print(" "+vis+": "+selfcal_plan[vis]['solint_settings'][selfcal_plan['solints'][iteration]]['interval']) + + # Set some cocal parameters. if selfcal_plan['solints'][iteration] in ["inf_EB_fb","inf_fb1"]: calculate_inf_EB_fb_anyways = True @@ -103,7 +107,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ calculate_inf_EB_fb_anyways = False preapply_targets_own_inf_EB = True # If there was no inf solint (e.g. because each source was observed only a single time, skip this as there are no gain tables to stick together. - if "inf" not in selfcal_plan['solints']: + if "inf" not in [selfcal_plan[vislist[0]]['solint_settings'][solint]['sub-name'] for solint in selfcal_plan['solints']]: iteration += 1 continue @@ -113,23 +117,41 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ # make sure the last phase-only selfcal table gets pre-applied # solints that use combinespw don't get pre-apply label by default for vis in vislist: - selfcal_plan[vis]['solint_settings'][selfcal_library['final_phase_solint']]['preapply_this_gaintable']=True + selfcal_plan[vis]['solint_settings'][selfcal_library[vis]['final_phase_solint']]['preapply_this_gaintable']=True - - if mode == "selfcal" and selfcal_plan['solint_snr'][selfcal_plan['solints'][iteration]] < minsnr_to_proceed and np.all([selfcal_plan[fid]['solint_snr_per_field'][selfcal_plan['solints'][iteration]] < minsnr_to_proceed for fid in selfcal_library['sub-fields']]): - print('*********** estimated SNR for solint='+selfcal_plan['solints'][iteration]+' too low, measured: '+str(selfcal_plan['solint_snr'][selfcal_plan['solints'][iteration]])+', Min SNR Required: '+str(minsnr_to_proceed)+' **************') + if mode == "selfcal": + remove_vis = [] + for vis in vislist: + if selfcal_plan[vis]['solint_snr'][selfcal_plan['solints'][iteration]] < minsnr_to_proceed and \ + np.all([selfcal_plan[fid][vis]['solint_snr_per_field'][selfcal_plan['solints'][iteration]] < minsnr_to_proceed for fid in \ + selfcal_library['sub-fields']]): + print('*********** estimated SNR for EB='+vis+' for solint='+selfcal_plan['solints'][iteration]+' too low, measured: '+\ + str(selfcal_plan[vis]['solint_snr'][selfcal_plan['solints'][iteration]])+', Min SNR Required: '+str(minsnr_to_proceed)+\ + ' **************') + remove_vis.append(vis) + + selfcal_library[vis]['Stop_Reason']='Estimated_SNR_too_low_for_solint '+selfcal_plan['solints'][iteration] + for fid in selfcal_library['sub-fields-to-selfcal']: + selfcal_library[fid][vis]['Stop_Reason']='Estimated_SNR_too_low_for_solint '+selfcal_plan['solints'][iteration] + + for rvis in remove_vis: + vislist.remove(rvis) + + if len(vislist) == 0: + print('*********** estimated SNR for solint='+selfcal_plan['solints'][iteration]+' too low for all EBs **************') if iteration > 1 and selfcal_plan['solmode'][iteration] !='ap' and do_amp_selfcal: # if a solution interval shorter than inf for phase-only SC has passed, attempt amplitude selfcal iteration=selfcal_plan['solmode'].index('ap') print('****************Attempting amplitude selfcal*************') continue - selfcal_library['Stop_Reason']='Estimated_SNR_too_low_for_solint '+selfcal_plan['solints'][iteration] - for fid in selfcal_library['sub-fields-to-selfcal']: - selfcal_library[fid]['Stop_Reason']='Estimated_SNR_too_low_for_solint '+selfcal_plan['solints'][iteration] break else: + selfcal_library['vislist-to-gaincal'] = vislist + for fid in selfcal_library['sub-fields']: + selfcal_library[fid]['vislist-to-gaincal'] = [vis for vis in selfcal_library['vislist-to-gaincal'] if vis in + selfcal_library[fid]['vislist']] + solint=selfcal_plan['solints'][iteration] - solint_interval=selfcal_plan['solint_interval'][iteration] if iteration == 0: print('Starting with solint: '+solint) else: @@ -144,7 +166,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ ## if selfcal_library['final_solint'] != 'None': prev_solint = selfcal_library['final_solint'] - prev_iteration = selfcal_library[vislist[0]][prev_solint]['iteration'] + prev_iteration = selfcal_library[[vis for vis in selfcal_library['vislist-to-gaincal'] if prev_solint in selfcal_library[vis]][0]][prev_solint]['iteration'] nterms_changed = (len(glob.glob(sani_target+'_'+band+'_'+prev_solint+'_'+str(prev_iteration)+"_post.model.tt*")) < selfcal_library['nterms']) @@ -168,7 +190,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ selfcal_library[vis][solint]={} selfcal_library[vis][solint]['clean_threshold'] = selfcal_library['nsigma'][iteration]*selfcal_library['RMS_NF_curr'] selfcal_library[vis][solint]['nfrms_multiplier'] = selfcal_library['RMS_NF_curr'] / selfcal_library['RMS_curr'] - for fid in np.intersect1d(selfcal_library['sub-fields-to-selfcal'],list(selfcal_library['sub-fields-fid_map'][vis].keys())): + for fid in np.intersect1d(selfcal_library[vis]['sub-fields-to-selfcal'],list(selfcal_library['sub-fields-fid_map'][vis].keys())): selfcal_library[fid][vis][solint]={} selfcal_library[fid][vis][solint]['clean_threshold'] = selfcal_library['nsigma'][iteration]*selfcal_library['RMS_NF_curr'] selfcal_library[fid][vis][solint]['nfrms_multiplier'] = selfcal_library['RMS_NF_curr'] / selfcal_library['RMS_curr'] @@ -185,13 +207,12 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ # Check that a mask was actually created, because if not the model will be empty and gaincal will do bad things and the # code will break. if not checkmask(sani_target+'_'+band+'_'+solint+'_'+str(iteration)+'.image.tt0') and mode !="cocal": - selfcal_library['Stop_Reason'] = 'Empty model for solint '+solint - for fid in selfcal_library['sub-fields-to-selfcal']: - selfcal_library[fid]['Stop_Reason'] = 'Empty model for solint '+solint for vis in vislist: + selfcal_library[vis]['Stop_Reason'] = 'Empty model for solint '+solint selfcal_library[vis][solint]['Pass'] = False selfcal_library[vis][solint]['Fail_Reason'] = 'Empty model for solint '+solint - for fid in np.intersect1d(selfcal_library['sub-fields-to-selfcal'],list(selfcal_library['sub-fields-fid_map'][vis].keys())): + for fid in np.intersect1d(selfcal_library[vis]['sub-fields-to-selfcal'],list(selfcal_library['sub-fields-fid_map'][vis].keys())): + selfcal_library[fid][vis]['Stop_Reason'] = 'Empty model for solint '+solint selfcal_library[fid][vis][solint]['Pass'] = False selfcal_library[fid][vis][solint]['Fail_Reason'] = 'Empty model for solint '+solint break # breakout of loop because the model is empty and gaincal will therefore fail @@ -224,25 +245,33 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ # Fields that don't have any mask in the primary beam should be removed from consideration, as their models are likely bad. if selfcal_library['obstype'] == 'mosaic': - selfcal_library['sub-fields-to-gaincal'] = evaluate_subfields_to_gaincal(selfcal_library, target, band, - solint, iteration, selfcal_plan['solmode'], selfcal_plan['solints'], selfcal_plan, minsnr_to_proceed, - allow_gain_interpolation=allow_gain_interpolation) + selfcal_library['sub-fields-to-gaincal'] = [] + for vis in vislist: + selfcal_library[vis]['sub-fields-to-gaincal'] = evaluate_subfields_to_gaincal(vis, selfcal_library, target, band, + solint, iteration, selfcal_plan['solmode'], selfcal_plan['solints'], selfcal_plan, minsnr_to_proceed, + allow_gain_interpolation=allow_gain_interpolation) + selfcal_library['sub-fields-to-gaincal'] = np.union1d(selfcal_library['sub-fields-to-gaincal'], selfcal_library[vis]['sub-fields-to-gaincal']).tolist() if solint != 'inf_EB' and not allow_gain_interpolation: selfcal_library['sub-fields-to-selfcal'] = selfcal_library['sub-fields-to-gaincal'] - print('Fields to gaincal: ',selfcal_library['sub-fields-to-gaincal']) + for vis in vislist: + selfcal_library[vis]['sub-fields-to-selfcal'] = selfcal_library['sub-fields-to-gaincal'] + + print('Fields to gaincal: ') + for vis in vislist: + print(f'{vis}:', selfcal_library[vis]['sub-fields-to-gaincal']) + if len(selfcal_library['sub-fields-to-gaincal']) == 0: - print('No fields to selfcal, exiting solution interval an selfcal for current target') - selfcal_library['Stop_Reason']='Missing_flux_in_all_sub-fields_for_solint '+selfcal_plan['solints'][iteration] - for fid in list(selfcal_library['sub-fields-fid_map'][vis].keys()): - selfcal_library[fid]['Stop_Reason']='Missing_flux_in_all_sub-fields_for_solint '+selfcal_plan['solints'][iteration] + print('No fields to selfcal, exiting solution interval an selfcal for current target') for vis in vislist: + selfcal_library[vis]['Stop_Reason']='Missing_flux_in_all_sub-fields_for_solint '+selfcal_plan['solints'][iteration] selfcal_library[vis][solint]['Pass'] = 'None' selfcal_library[vis][solint]['Fail_Reason'] = 'Missing_flux_in_all_sub-fields_for_solint '+solint #for fid in np.intersect1d(selfcal_library['sub-fields-to-selfcal'],list(selfcal_library['sub-fields-fid_map'][vis].keys())): for fid in list(selfcal_library['sub-fields-fid_map'][vis].keys()): if solint in selfcal_library[fid][vis].keys(): + selfcal_library[fid][vis]['Stop_Reason']='Missing_flux_in_all_sub-fields_for_solint '+selfcal_plan['solints'][iteration] selfcal_library[fid][vis][solint]['Pass'] = 'None' selfcal_library[fid][vis][solint]['Fail_Reason'] = 'Missing_flux_in_all_sub-fields_for_solint '+solint @@ -251,17 +280,21 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ else: selfcal_library['sub-fields-to-gaincal'] = selfcal_library['sub-fields-to-selfcal'] + for vis in vislist: + selfcal_library[vis]['sub-fields-to-gaincal'] = selfcal_library[vis]['sub-fields-to-selfcal'] - + print('selfcal_library["sub-fields-to-selfcal"] = ', selfcal_library['sub-fields-to-selfcal']) + for vis in selfcal_library['vislist']: + print('selfcal_library[vis]["sub-fields-to-selfcal"] = ', selfcal_library[vis]['sub-fields-to-selfcal']) # Calculate the complex gains for vis in vislist: - if np.intersect1d(selfcal_library['sub-fields-to-gaincal'],\ + if np.intersect1d(selfcal_library[vis]['sub-fields-to-gaincal'],\ list(selfcal_library['sub-fields-fid_map'][vis].keys())).size == 0: continue - gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, solint_interval, selfcal_plan['applycal_mode'][iteration], \ + gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, selfcal_plan['applycal_mode'][iteration], \ iteration, gaincal_minsnr, \ gaincal_unflag_minsnr=gaincal_unflag_minsnr, minsnr_to_proceed=minsnr_to_proceed, rerank_refants=rerank_refants, \ unflag_only_lbants=unflag_only_lbants, unflag_only_lbants_onlyap=unflag_only_lbants_onlyap, calonly_max_flagged=calonly_max_flagged, @@ -273,16 +306,23 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ # With gaincal done and bad fields removed from gain tables if necessary, check whether any fields should no longer be # selfcal'd because they have too much interpolation. if selfcal_library['obstype'] == 'mosaic': - selfcal_library['sub-fields-to-selfcal'] = evaluate_subfields_after_gaincal(selfcal_library, target, band, - solint, iteration, selfcal_plan['solmode'], allow_gain_interpolation=allow_gain_interpolation) + selfcal_library['sub-fields-to-selfcal'] = [] + for vis in vislist: + selfcal_library[vis]['sub-fields-to-selfcal'] = evaluate_subfields_after_gaincal(vis, selfcal_library, selfcal_plan, target, band, + solint, iteration, selfcal_plan['solmode'], allow_gain_interpolation=allow_gain_interpolation) + selfcal_library['sub-fields-to-selfcal'] = np.union1d(selfcal_library['sub-fields-to-selfcal'], selfcal_library[vis]['sub-fields-to-selfcal']).astype(int).tolist() + + print('selfcal_library["sub-fields-to-selfcal"] = ', selfcal_library['sub-fields-to-selfcal']) + for vis in selfcal_library['vislist']: + print('selfcal_library[vis]["sub-fields-to-selfcal"] = ', selfcal_library[vis]['sub-fields-to-selfcal']) ## ## Apply gain solutions per MS, target, solint, and band ## for vis in vislist: applycal_wrapper(vis, target, band, solint, selfcal_library, - current=lambda f: f in selfcal_library['sub-fields-to-selfcal'], - final=lambda f: f not in selfcal_library['sub-fields-to-selfcal'] and selfcal_library[f]['SC_success'], + current=lambda f: f in selfcal_library[vis]['sub-fields-to-selfcal'], + final=lambda f: f not in selfcal_library[vis]['sub-fields-to-selfcal'] and selfcal_library[f]['SC_success'], restore_flags='fb_selfcal_starting_flags_'+sani_target+'_'+band if mode == "cocal" else None) ## Create post self-cal image using the model as a startmodel to evaluate how much selfcal helped @@ -312,6 +352,9 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ mosaic_SNR, mosaic_RMS, mosaic_SNR_NF, mosaic_RMS_NF = {}, {}, {}, {} post_mosaic_SNR, post_mosaic_RMS, post_mosaic_SNR_NF, post_mosaic_RMS_NF = {}, {}, {}, {} + print('selfcal_library["sub-fields-to-selfcal"] = ', selfcal_library['sub-fields-to-selfcal']) + for vis in vislist: + print('selfcal_library[vis]["sub-fields-to-selfcal"] = ', selfcal_library[vis]['sub-fields-to-selfcal']) for fid in selfcal_library['sub-fields-to-selfcal']: if selfcal_library['obstype'] == 'mosaic': imagename = sani_target+'_field_'+str(fid)+'_'+band+'_'+solint+'_'+str(iteration) @@ -374,9 +417,9 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ (post_mosaic_RMS_NF[fid] - mosaic_RMS_NF[fid])/mosaic_RMS_NF[fid] < 1.05) or \ (((post_mosaic_RMS[fid] - mosaic_RMS[fid])/mosaic_RMS[fid] > 1.05 or \ (post_mosaic_RMS_NF[fid] - mosaic_RMS_NF[fid])/mosaic_RMS_NF[fid] > 1.05) and \ - selfcal_plan[fid]['solint_snr_per_field'][solint] > 5) + np.all([selfcal_plan[fid][vis]['solint_snr_per_field'][solint] > 5 for vis in selfcal_library['vislist-to-gaincal']])) - if 'inf_EB' in solint or 'delay' in solint: + if 'inf_EB' in solint or 'd' in solint: # If any of the fields succeed in the "strict" sense, then allow for minor reductions in the evaluation quantity in other # fields because there's a good chance that those are just noise being pushed around. field_by_field_success = numpy.logical_and(numpy.logical_and(loose_field_by_field_success, beam_field_by_field_success), \ @@ -392,10 +435,10 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ print('****************Not all fields were successful, so re-applying and re-making _post image*************') for vis in vislist: applycal_wrapper(vis, target, band, solint, selfcal_library, - current=lambda f: f in field_by_field_success_dict and field_by_field_success_dict[f], - final=lambda f: (f not in field_by_field_success_dict or not field_by_field_success_dict[f]) and + current=lambda f: f in field_by_field_success_dict and field_by_field_success_dict[f] and f in selfcal_library[vis]['sub-fields-to-selfcal'], + final=lambda f: (f not in field_by_field_success_dict or not field_by_field_success_dict[f] or f not in selfcal_library[vis]['sub-fields-to-selfcal']) and selfcal_library[f]['SC_success'], - clear=lambda f: (f not in field_by_field_success_dict or not field_by_field_success_dict[f]) and + clear=lambda f: (f not in field_by_field_success_dict or not field_by_field_success_dict[f] or f not in selfcal_library[vis]['sub-fields-to-selfcal']) and not selfcal_library[f]['SC_success'], restore_flags='selfcal_starting_flags_'+sani_target+'_'+band) @@ -454,23 +497,24 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ #run a pre-check as to whether a marginal inf_EB result will go on to attempt inf, if not we will fail a marginal inf_EB marginal_inf_EB_will_attempt_next_solint=False - if (solint =='inf_EB' or 'delay' in solint) and ((((post_SNR-SNR)/SNR > -0.02) and ((post_SNR-SNR)/SNR < 0.00)) or (((post_SNR_NF - SNR_NF)/SNR_NF > -0.02) and ((post_SNR_NF - SNR_NF)/SNR_NF < 0.00))) and (delta_beamarea < delta_beam_thresh): - if selfcal_plan['solint_snr'][selfcal_plan['solints'][iteration+1]] < minsnr_to_proceed and np.all([selfcal_plan[fid]['solint_snr_per_field'][selfcal_plan['solints'][iteration+1]] < minsnr_to_proceed for fid in selfcal_library['sub-fields']]): + if (solint =='inf_EB' or 'd' in solint) and ((((post_SNR-SNR)/SNR > -0.02) and ((post_SNR-SNR)/SNR < 0.00)) or (((post_SNR_NF - SNR_NF)/SNR_NF > -0.02) and ((post_SNR_NF - SNR_NF)/SNR_NF < 0.00))) and (delta_beamarea < delta_beam_thresh): + if np.all([selfcal_plan[vis]['solint_snr'][selfcal_plan['solints'][iteration+1]] < minsnr_to_proceed for vis in selfcal_library['vislist-to-gaincal']]) and \ + np.all([np.all([selfcal_plan[fid][vis]['solint_snr_per_field'][selfcal_plan['solints'][iteration+1]] < minsnr_to_proceed for fid in selfcal_library['sub-fields']]) for vis in selfcal_library['vislist-to-gaincal']]): marginal_inf_EB_will_attempt_next_solint = False else: marginal_inf_EB_will_attempt_next_solint = True RMS_change_acceptable= False if mode != 'cocal': RMS_change_acceptable = (post_RMS/RMS < 1.05 and post_RMS_NF/RMS_NF < 1.05) or \ - ((post_RMS/RMS > 1.05 or post_RMS_NF/RMS_NF > 1.05) and selfcal_plan['solint_snr'][solint] > 5) + ((post_RMS/RMS > 1.05 or post_RMS_NF/RMS_NF > 1.05) and np.any([selfcal_plan[vis]['solint_snr'][solint] > 5 for vis in selfcal_library['vislist-to-gaincal']])) else: RMS_change_acceptable = (post_RMS/RMS < 1.05 and post_RMS_NF/RMS_NF < 1.05) - if (((post_SNR >= SNR) and (post_SNR_NF >= SNR_NF) and (delta_beamarea < delta_beam_thresh)) or ((('inf_EB' in solint) or 'delay' in solint) and marginal_inf_EB_will_attempt_next_solint and ((post_SNR-SNR)/SNR > -0.02) and ((post_SNR_NF - SNR_NF)/SNR_NF > -0.02) and (delta_beamarea < delta_beam_thresh))) and np.any(field_by_field_success) and RMS_change_acceptable: + if (((post_SNR >= SNR) and (post_SNR_NF >= SNR_NF) and (delta_beamarea < delta_beam_thresh)) or ((('inf_EB' in solint) or 'd' in solint) and marginal_inf_EB_will_attempt_next_solint and ((post_SNR-SNR)/SNR > -0.02) and ((post_SNR_NF - SNR_NF)/SNR_NF > -0.02) and (delta_beamarea < delta_beam_thresh))) and np.any(field_by_field_success) and RMS_change_acceptable: if do_fallback_combinespw: for vis in vislist: - if 'delay' in solint: + if 'd' in solint: selfcal_plan[vis]['solint_settings'][solint]['final_mode']='per_bb' else: selfcal_plan[vis]['solint_settings'][solint]['final_mode']='combinespw' @@ -493,13 +537,13 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ continue selfcal_library['SC_success']=True - selfcal_library['Stop_Reason']='None' #keep track of whether inf_EB had a S/N decrease - if (solint =='inf_EB' or 'delay' in solint) and (((post_SNR-SNR)/SNR < 0.0) or ((post_SNR_NF - SNR_NF)/SNR_NF < 0.0)): + if (solint =='inf_EB' or 'd' in solint) and (((post_SNR-SNR)/SNR < 0.0) or ((post_SNR_NF - SNR_NF)/SNR_NF < 0.0)): selfcal_library['inf_EB_SNR_decrease']=True - elif (solint =='inf_EB' or 'delay' in solint) and (((post_SNR-SNR)/SNR >= 0.0) and ((post_SNR_NF - SNR_NF)/SNR_NF >= 0.0)): + elif (solint =='inf_EB' or 'd' in solint) and (((post_SNR-SNR)/SNR >= 0.0) and ((post_SNR_NF - SNR_NF)/SNR_NF >= 0.0)): selfcal_library['inf_EB_SNR_decrease']=False for vis in vislist: + selfcal_library[vis]['Stop_Reason']='None' selfcal_library[vis]['gaintable_final']=selfcal_library[vis][solint]['gaintable'] selfcal_library[vis]['spwmap_final']=selfcal_library[vis][solint]['spwmap'].copy() selfcal_library[vis]['applycal_mode_final']=selfcal_library[vis][solint]['applycal_mode'] @@ -507,8 +551,9 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ selfcal_library[vis]['gaincal_combine_final']=selfcal_library[vis][solint]['gaincal_combine'] selfcal_library[vis][solint]['Pass']=True selfcal_library[vis][solint]['Fail_Reason']='None' - if selfcal_plan['solmode'][iteration]=='p': - selfcal_library['final_phase_solint']=solint + if selfcal_plan['solmode'][iteration]=='p': + selfcal_library[vis]['final_phase_solint']=solint + selfcal_library[vis]['final_solint']=solint selfcal_library['final_solint']=solint selfcal_library['final_solint_mode']=selfcal_plan['solmode'][iteration] @@ -517,13 +562,19 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ for ind, fid in enumerate(selfcal_library['sub-fields-to-selfcal']): if field_by_field_success[ind]: selfcal_library[fid]['SC_success']=True - selfcal_library[fid]['Stop_Reason']='None' if (solint =='inf_EB') and not strict_field_by_field_success[ind]: selfcal_library[fid]['inf_EB_SNR_decrease']=True elif (solint =='inf_EB') and strict_field_by_field_success[ind]: selfcal_library[fid]['inf_EB_SNR_decrease']=False - for vis in selfcal_library[fid]['vislist']: + for vis in selfcal_library[fid]['vislist-to-gaincal']: + if fid not in selfcal_library[vis]['sub-fields-to-selfcal']: + continue + # I think the below might not be correct for mosaics - it would set the gaintable even if fid is not in [vis]['sub-fields-to-selfcal']. + # I think this needs: + # if fid not in selfcal_library[vis]['sub-fields-to-selfcal']: + # continue + selfcal_library[fid][vis]['Stop_Reason']='None' selfcal_library[fid][vis]['gaintable_final']=selfcal_library[fid][vis][solint]['gaintable'] selfcal_library[fid][vis]['spwmap_final']=selfcal_library[fid][vis][solint]['spwmap'].copy() selfcal_library[fid][vis]['applycal_mode_final']=selfcal_library[fid][vis][solint]['applycal_mode'] @@ -531,8 +582,9 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ selfcal_library[fid][vis]['gaincal_combine_final']=selfcal_library[fid][vis][solint]['gaincal_combine'] selfcal_library[fid][vis][solint]['Pass']=True selfcal_library[fid][vis][solint]['Fail_Reason']='None' - if selfcal_plan['solmode'][iteration]=='p': - selfcal_library[fid]['final_phase_solint']=solint + if selfcal_plan['solmode'][iteration]=='p': + selfcal_library[fid][vis]['final_phase_solint']=solint + selfcal_library[fid][vis]['final_solint']=solint selfcal_library[fid]['final_solint']=solint selfcal_library[fid]['final_solint_mode']=selfcal_plan['solmode'][iteration] selfcal_library[fid]['iteration']=iteration @@ -600,7 +652,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ selfcal_library[vis][solint]['Pass']=False for fid in selfcal_library['sub-fields-to-selfcal']: - for vis in selfcal_library[fid]['vislist']: + for vis in selfcal_library[fid]['vislist-to-gaincal']: selfcal_library[fid][vis][solint]['Pass']=False repeat_solint = False @@ -613,7 +665,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ ## if S/N worsens, and/or beam area increases reject current solutions and reapply previous (or revert to origional data) ## - if not selfcal_library[vislist[0]][solint]['Pass'] or (solint == 'inf_EB' and selfcal_library['inf_EB_SNR_decrease']): + if not selfcal_library[selfcal_library['vislist-to-gaincal'][0]][solint]['Pass'] or (solint == 'inf_EB' and selfcal_library['inf_EB_SNR_decrease']): reason='' if (post_SNR <= SNR): reason=reason+' S/N decrease' @@ -627,11 +679,11 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ reason=reason+'Beam change beyond '+str(delta_beam_thresh) if mode != 'cocal': - if (post_RMS/RMS > 1.05 and selfcal_plan['solint_snr'][solint] <= 5): + if (post_RMS/RMS > 1.05 and np.all([selfcal_plan[vis]['solint_snr'][solint] <= 5 for vis in selfcal_library['vislist-to-gaincal']])): if reason != '': reason=reason+'; ' reason=reason+'RMS increase beyond 5%' - if (post_RMS_NF/RMS_NF > 1.05 and selfcal_plan['solint_snr'][solint] <= 5): + if (post_RMS_NF/RMS_NF > 1.05 and np.all([selfcal_plan[vis]['solint_snr'][solint] <= 5 for vis in selfcal_library['vislist-to-gaincal']])): if reason != '': reason=reason+'; ' reason=reason+'NF RMS increase beyond 5%' @@ -649,14 +701,14 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ if reason != '': reason=reason+'; ' reason=reason+'All sub-fields failed' - selfcal_library['Stop_Reason']=reason for vis in vislist: + selfcal_library[vis]['Stop_Reason']=reason #selfcal_library[vis][solint]['Pass']=False selfcal_library[vis][solint]['Fail_Reason']=reason mosaic_reason = {} for fid in selfcal_library['sub-fields-to-selfcal']: - if not selfcal_library[fid][selfcal_library[fid]['vislist'][0]][solint]['Pass'] or \ + if not selfcal_library[fid][selfcal_library[fid]['vislist-to-gaincal'][0]][solint]['Pass'] or \ (solint == "inf_EB" and selfcal_library[fid]['inf_EB_SNR_decrease']): mosaic_reason[fid]='' if (post_mosaic_SNR[fid] <= mosaic_SNR[fid]): @@ -669,24 +721,24 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ if mosaic_reason[fid] !='': mosaic_reason[fid]=mosaic_reason[fid]+'; ' mosaic_reason[fid]=mosaic_reason[fid]+'Beam change beyond '+str(delta_beam_thresh) - if (post_RMS/RMS > 1.05 and selfcal_plan['solint_snr'][solint] <= 5): + if (post_RMS/RMS > 1.05 and np.all([selfcal_plan[vis]['solint_snr'][solint] <= 5 for vis in selfcal_library[fid]['vislist-to-gaincal']])): if mosaic_reason[fid] != '': mosaic_reason[fid]=mosaic_reason[fid]+'; ' mosaic_reason[fid]=mosaic_reason[fid]+'RMS increase beyond 5%' - if (post_RMS_NF/RMS_NF > 1.05 and selfcal_plan['solint_snr'][solint] <= 5): + if (post_RMS_NF/RMS_NF > 1.05 and np.all([selfcal_plan[vis]['solint_snr'][solint] <= 5 for vis in selfcal_library[fid]['vislist-to-gaincal']])): if mosaic_reason[fid] != '': mosaic_reason[fid]=mosaic_reason[fid]+'; ' mosaic_reason[fid]=mosaic_reason[fid]+'NF RMS increase beyond 5%' if mosaic_reason[fid] == '': mosaic_reason[fid] = "Global selfcal failed" - selfcal_library[fid]['Stop_Reason']=mosaic_reason[fid] - for vis in selfcal_library[fid]['vislist']: + for vis in selfcal_library[fid]['vislist-to-gaincal']: + selfcal_library[fid][vis]['Stop_Reason']=mosaic_reason[fid] #selfcal_library[fid][vis][solint]['Pass']=False selfcal_library[fid][vis][solint]['Fail_Reason']=mosaic_reason[fid] # If any of the fields failed self-calibration, we need to re-apply calibrations for all fields because we need to revert flagging back # to the starting point. - if np.any([selfcal_library[fid][selfcal_library[fid]['vislist'][0]][solint]['Pass'] == False for fid in \ + if np.any([selfcal_library[fid][selfcal_library[fid]['vislist-to-gaincal'][0]][solint]['Pass'] == False for fid in \ selfcal_library['sub-fields-to-selfcal']]) or len(selfcal_library['sub-fields-to-selfcal']) < \ len(selfcal_library['sub-fields']): print('****************Selfcal failed for some sub-fields:*************') @@ -714,13 +766,17 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ # because the solint failed, because for mosaics we can keep trying the field as we clean deeper. If we set to False now, that wont happen. for fid in np.intersect1d(selfcal_library['sub-fields'],list(selfcal_library['sub-fields-fid_map'][vis].keys())): if (selfcal_library['final_solint'] == 'inf_EB' and selfcal_library['inf_EB_SNR_decrease']) or \ - ((not selfcal_library[vislist[0]][solint]['Pass'] or solint == 'int') and \ + (np.any([not selfcal_library[vis][solint]['Pass'] or selfcal_plan[vis]['solint_settings'][solint]['sub-name'] == 'int' for vis in selfcal_library[fid]['vislist-to-gaincal']]) and \ (selfcal_library[fid]['final_solint'] == 'inf_EB' and selfcal_library[fid]['inf_EB_SNR_decrease'])): - selfcal_library[fid]['SC_success']=False - selfcal_library[fid]['final_solint']='None' - for vis in selfcal_library[fid]['vislist']: - selfcal_library[fid][vis]['inf_EB']['Pass']=False # remove the success from inf_EB - selfcal_library[fid][vis]['inf_EB']['Fail_Reason']+=' with no successful solints later' # remove the success from inf_EB + + for vis in selfcal_library[fid]['vislist-to-gaincal']: + if not selfcal_library[vis][solint]['Pass'] or selfcal_plan[vis]['solint_settings'][solint]['sub-name'] == 'int': + selfcal_library[fid][vis]['inf_EB']['Pass']=False # remove the success from inf_EB + selfcal_library[fid][vis]['inf_EB']['Fail_Reason']+=' with no successful solints later' # remove the success from inf_EB + + if np.all([not selfcal_library[fid][vis]['inf_EB']['Pass'] for vis in selfcal_library[fid]['vislist']]): + selfcal_library[fid]['SC_success']=False + selfcal_library[fid]['final_solint']='None' for vis in vislist: applycal_wrapper(vis, target, band, solint, selfcal_library, @@ -739,25 +795,31 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ # If any of the sub-fields passed, and the whole mosaic passed, then we can move on to the next solint, otherwise we have to back out. if selfcal_library[vislist[0]][solint]['Pass'] == True and \ - np.any([selfcal_library[fid][selfcal_library[fid]['vislist'][0]][solint]['Pass'] == True for fid in \ + np.any([selfcal_library[fid][selfcal_library[fid]['vislist-to-gaincal'][0]][solint]['Pass'] == True for fid in \ selfcal_library['sub-fields-to-selfcal']]): if mode == "selfcal" and (iteration < len(selfcal_plan['solints'])-1) and (selfcal_library[vis][solint]['SNR_post'] > \ selfcal_library['SNR_orig']): #(iteration == 0) and print('Updating solint = '+selfcal_plan['solints'][iteration+1]+' SNR') - print('Was: ',selfcal_plan['solint_snr'][selfcal_plan['solints'][iteration+1]]) - get_SNR_self_update(selfcal_library,selfcal_plan,n_ants,solint,selfcal_plan['solints'][iteration+1],selfcal_plan['integration_time'], - selfcal_plan['solint_snr']) - print('Now: ',selfcal_plan['solint_snr'][selfcal_plan['solints'][iteration+1]]) + for vis in vislist: + if selfcal_plan['solints'][iteration+1] not in selfcal_plan[vis]['solint_snr']: + continue + print(vis+' was: ',selfcal_plan[vis]['solint_snr'][selfcal_plan['solints'][iteration+1]]) + get_SNR_self_update(vis, selfcal_library,selfcal_plan,n_ants,solint,selfcal_plan['solints'][iteration+1],selfcal_plan[vis]['integration_time'], + selfcal_plan[vis]['solint_snr']) + print(vis+' now: ',selfcal_plan[vis]['solint_snr'][selfcal_plan['solints'][iteration+1]]) for fid in selfcal_library['sub-fields-to-selfcal']: - print('Field '+str(fid)+' Was: ',selfcal_plan[fid]['solint_snr_per_field'][selfcal_plan['solints'][iteration+1]]) - get_SNR_self_update(selfcal_library[fid],selfcal_plan,n_ants,solint,selfcal_plan['solints'][iteration+1], - selfcal_plan['integration_time'],selfcal_plan[fid]['solint_snr_per_field']) - print('Field '+str(fid)+' Now: ',selfcal_plan[fid]['solint_snr_per_field'][selfcal_plan['solints'][iteration+1]]) + for vis in selfcal_library[fid]['vislist']: + if selfcal_plan['solints'][iteration+1] not in selfcal_plan[fid][vis]['solint_snr_per_field']: + continue + print('Field '+str(fid)+' '+vis+' was: ',selfcal_plan[fid][vis]['solint_snr_per_field'][selfcal_plan['solints'][iteration+1]]) + get_SNR_self_update(vis, selfcal_library[fid],selfcal_plan,n_ants,solint,selfcal_plan['solints'][iteration+1], + selfcal_plan[vis]['integration_time'],selfcal_plan[fid][vis]['solint_snr_per_field']) + print('Field '+str(fid)+' '+vis+' now: ',selfcal_plan[fid][vis]['solint_snr_per_field'][selfcal_plan['solints'][iteration+1]]) # If not all fields succeed for inf_EB or scan_inf/inf, depending on mosaic or single field, then don't go on to amplitude selfcal, # even if *some* fields succeeded. - if iteration <= 1 and ((not np.all([selfcal_library[fid][selfcal_library[fid]['vislist'][0]][solint]['Pass'] == True for fid in \ + if iteration <= 1 and ((not np.all([selfcal_library[fid][selfcal_library[fid]['vislist-to-gaincal'][0]][solint]['Pass'] == True for fid in \ selfcal_library['sub-fields-to-selfcal']])) or len(selfcal_library['sub-fields-to-selfcal']) < \ len(selfcal_library['sub-fields'])) and do_amp_selfcal: print("***** NOTE: Amplitude self-calibration turned off because not all fields succeeded at non-inf_EB phase self-calibration") @@ -779,6 +841,12 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ if mode == "selfcal" and iteration > 1 and selfcal_plan['solmode'][iteration] !='ap' and do_amp_selfcal: # if a solution interval shorter than inf for phase-only SC has passed, attempt amplitude selfcal iteration=selfcal_plan['solmode'].index('ap') selfcal_library['sub-fields-to-selfcal'] = selfcal_library['sub-fields'] + for vis in vislist: + selfcal_library[vis]['sub-fields-to-selfcal'] = selfcal_library['sub-fields'].copy() + for fid in selfcal_library[vis]['sub-fields-to-selfcal']: + if not fid in selfcal_library['sub-fields-fid_map'][vis]: + # If a sub-field is not in an EB, it shouldn't be considered for selfcal for that EB + selfcal_library[vis]['sub-fields-to-selfcal'].pop(fid) print('****************Selfcal halted for phase, attempting amplitude*************') continue elif mode == "cocal" and "inf_fb" in solint: @@ -790,14 +858,20 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ break # breakout of loops of successive solints since solutions are getting worse # Finally, update the list of fields to be self-calibrated now that we don't need to know the list at the beginning of this solint. - new_fields_to_selfcal = [] coarsest_solint=selfcal_plan['solints'][0] - for fid in selfcal_library['sub-fields']: - if mode == "cocal": - if ("inf_EB" in selfcal_library[fid]['vislist'][0] and selfcal_library[fid][selfcal_library[fid]['vislist'][0]][coarsest_solint]["Pass"]) or selfcal_library[fid][selfcal_library[fid]['vislist'][0]][coarsest_solint+"_fb"]["Pass"]: - new_fields_to_selfcal.append(fid) - else: - if selfcal_library[fid][selfcal_library[fid]['vislist'][0]][coarsest_solint]["Pass"]: - new_fields_to_selfcal.append(fid) + selfcal_library['sub-fields-to-selfcal'] = [] + for vis in vislist: + new_fields_to_selfcal = [] + for fid in selfcal_library['sub-fields']: + if vis not in selfcal_library[fid]['vislist']: + continue + + if mode == "cocal": + if ("inf_EB" in selfcal_library[fid]['vislist'][0] and selfcal_library[fid][vis][coarsest_solint]["Pass"]) or selfcal_library[fid][selfcal_library[fid]['vislist'][0]][coarsest_solint+"_fb"]["Pass"]: + new_fields_to_selfcal.append(fid) + else: + if selfcal_library[fid][vis][coarsest_solint]["Pass"]: + new_fields_to_selfcal.append(fid) - selfcal_library['sub-fields-to-selfcal'] = new_fields_to_selfcal + selfcal_library[vis]['sub-fields-to-selfcal'] = new_fields_to_selfcal + selfcal_library['sub-fields-to-selfcal'] = np.union1d(selfcal_library['sub-fields-to-selfcal'], selfcal_library[vis]['sub-fields-to-selfcal']).tolist() diff --git a/auto_selfcal/selfcal_helpers.py b/auto_selfcal/selfcal_helpers.py index 4dc9d20..1c40721 100644 --- a/auto_selfcal/selfcal_helpers.py +++ b/auto_selfcal/selfcal_helpers.py @@ -1787,43 +1787,18 @@ def get_SNR_self(selfcal_library,selfcal_plan,n_ant,inf_EB_gaincal_combine,inf_E minsolint_spw=100 for target in selfcal_library: for band in selfcal_library[target].keys(): - selfcal_plan[target][band]['solint_snr'], selfcal_plan[target][band]['solint_snr_per_spw'], selfcal_plan[target][band]['solint_snr_per_bb'] = \ - get_SNR_self_individual(selfcal_library[target][band]['vislist'], selfcal_library[target][band], n_ant, selfcal_plan[target][band]['solints'], - selfcal_plan[target][band]['solint_interval'], - selfcal_plan[target][band]['integration_time'], inf_EB_gaincal_combine, inf_EB_gaintype) - - print('Estimated SNR per solint:') - print(target,band) - for solint in selfcal_plan[target][band]['solints']: - if solint == 'inf_EB': - print('{}: {:0.2f}'.format(solint,selfcal_plan[target][band]['solint_snr'][solint])) - ''' - for spw in solint_snr_per_spw[target][band][solint].keys(): - print('{}: spw: {}: {:0.2f}, BW: {} GHz'.format(solint,spw,solint_snr_per_spw[target][band][solint][spw],selfcal_library[target][band]['per_spw_stats'][str(spw)]['effective_bandwidth'])) - if solint_snr_per_spw[target][band][solint][spw] < minsolint_spw: - minsolint_spw=solint_snr_per_spw[target][band][solint][spw] - if minsolint_spw < 3.5 and minsolint_spw > 2.5 and inf_EB_override==False: # if below 3.5 but above 2.5 switch to gaintype T, but leave combine=scan - print('Switching Gaintype to T for: '+target) - inf_EB_gaintype_dict[target][band]='T' - elif minsolint_spw < 2.5 and inf_EB_override==False: - print('Switching Gaincal combine to spw,scan for: '+target) - inf_EB_gaincal_combine_dict[target][band]='scan,spw' # if below 2.5 switch to combine=spw to avoid losing spws - ''' - else: - print('{}: {:0.2f}'.format(solint,selfcal_plan[target][band]['solint_snr'][solint])) - - for fid in selfcal_library[target][band]['sub-fields']: - selfcal_plan[target][band][fid] = {} - selfcal_plan[target][band][fid]['solint_snr_per_field'], selfcal_plan[target][band][fid]['solint_snr_per_field_per_spw'], selfcal_plan[target][band][fid]['solint_snr_per_field_per_bb'] = \ - get_SNR_self_individual(selfcal_library[target][band]['vislist'], selfcal_library[target][band][fid], n_ant, - selfcal_plan[target][band]['solints'], selfcal_plan[target][band]['solint_interval'], selfcal_plan[target][band]['integration_time'], inf_EB_gaincal_combine, - inf_EB_gaintype) + for vis in selfcal_library[target][band]['vislist']: + solints_per_vis = [solint for solint in selfcal_plan[target][band]['solints'] if solint in selfcal_plan[target][band][vis]['solint_settings']] + selfcal_plan[target][band][vis]['solint_snr'], selfcal_plan[target][band][vis]['solint_snr_per_spw'], selfcal_plan[target][band][vis]['solint_snr_per_bb'] = \ + get_SNR_self_individual([vis], selfcal_library[target][band], n_ant, solints_per_vis, + selfcal_plan[target][band][vis]['solint_settings'],selfcal_plan[target][band][vis]['integration_time'], inf_EB_gaincal_combine, inf_EB_gaintype) + print('Estimated SNR per solint:') - print(target,band,"field "+str(fid)) - for solint in selfcal_plan[target][band]['solints']: - if solint == 'inf_EB': - print('{}: {:0.2f}'.format(solint,selfcal_plan[target][band][fid]['solint_snr_per_field'][solint])) + print(target,band,vis) + for solint in solints_per_vis: + if selfcal_plan[target][band][vis]['solint_settings'][solint]['interval'] == 'inf_EB': + print('{}: {:0.2f}'.format(solint,selfcal_plan[target][band][vis]['solint_snr'][solint])) ''' for spw in solint_snr_per_spw[target][band][solint].keys(): print('{}: spw: {}: {:0.2f}, BW: {} GHz'.format(solint,spw,solint_snr_per_spw[target][band][solint][spw],selfcal_library[target][band]['per_spw_stats'][str(spw)]['effective_bandwidth'])) @@ -1837,11 +1812,43 @@ def get_SNR_self(selfcal_library,selfcal_plan,n_ant,inf_EB_gaincal_combine,inf_E inf_EB_gaincal_combine_dict[target][band]='scan,spw' # if below 2.5 switch to combine=spw to avoid losing spws ''' else: - print('{}: {:0.2f}'.format(solint,selfcal_plan[target][band][fid]['solint_snr_per_field'][solint])) + print('{}: {:0.2f}'.format(solint,selfcal_plan[target][band][vis]['solint_snr'][solint])) + + for fid in selfcal_library[target][band]['sub-fields']: + selfcal_plan[target][band][fid] = {} + for vis in selfcal_library[target][band][fid]['vislist']: + print(target, band, vis, fid, selfcal_library[target][band][fid]['vislist']) + selfcal_plan[target][band][fid][vis] = {} + solints_per_vis = [solint for solint in selfcal_plan[target][band]['solints'] if solint in selfcal_plan[target][band][vis]['solint_settings']] + selfcal_plan[target][band][fid][vis]['solint_snr_per_field'], selfcal_plan[target][band][fid][vis]['solint_snr_per_field_per_spw'], \ + selfcal_plan[target][band][fid][vis]['solint_snr_per_field_per_bb'] = \ + get_SNR_self_individual([vis], selfcal_library[target][band][fid], n_ant, + solints_per_vis, selfcal_plan[target][band][vis]['solint_settings'], \ + selfcal_plan[target][band][vis]['integration_time'], inf_EB_gaincal_combine, inf_EB_gaintype) + + print('Estimated SNR per solint:') + print(target,band,"field "+str(fid),vis) + for solint in solints_per_vis: + if selfcal_plan[target][band][vis]['solint_settings'][solint]['interval'] == 'inf_EB': + print('{}: {:0.2f}'.format(solint,selfcal_plan[target][band][fid][vis]['solint_snr_per_field'][solint])) + ''' + for spw in solint_snr_per_spw[target][band][solint].keys(): + print('{}: spw: {}: {:0.2f}, BW: {} GHz'.format(solint,spw,solint_snr_per_spw[target][band][solint][spw],selfcal_library[target][band]['per_spw_stats'][str(spw)]['effective_bandwidth'])) + if solint_snr_per_spw[target][band][solint][spw] < minsolint_spw: + minsolint_spw=solint_snr_per_spw[target][band][solint][spw] + if minsolint_spw < 3.5 and minsolint_spw > 2.5 and inf_EB_override==False: # if below 3.5 but above 2.5 switch to gaintype T, but leave combine=scan + print('Switching Gaintype to T for: '+target) + inf_EB_gaintype_dict[target][band]='T' + elif minsolint_spw < 2.5 and inf_EB_override==False: + print('Switching Gaincal combine to spw,scan for: '+target) + inf_EB_gaincal_combine_dict[target][band]='scan,spw' # if below 2.5 switch to combine=spw to avoid losing spws + ''' + else: + print('{}: {:0.2f}'.format(solint,selfcal_plan[target][band][fid][vis]['solint_snr_per_field'][solint])) #return solint_snr, solint_snr_per_spw, solint_snr_per_field, solint_snr_per_field_per_spw -def get_SNR_self_individual(vislist,selfcal_library,n_ant,solints,solints_interval,integration_time,inf_EB_gaincal_combine,inf_EB_gaintype): +def get_SNR_self_individual(vislist,selfcal_library,n_ant,solints,solint_settings,integration_time,inf_EB_gaincal_combine,inf_EB_gaintype): if inf_EB_gaintype=='G': polscale=2.0 else: @@ -1857,106 +1864,125 @@ def get_SNR_self_individual(vislist,selfcal_library,n_ant,solints,solints_interv solint_snr_per_spw[solint]={} solint_snr_per_bb[solint]={} if solint == 'inf_EB' or solint == 'inf_EB_delay': - SNR_self_EB=np.zeros(len(selfcal_library['vislist'])) + SNR_self_EB=np.zeros(len(vislist)) SNR_self_EB_spw={} SNR_self_EB_bb={} - for i in range(len(selfcal_library['vislist'])): - if solints_interval[s] != 'inf': - solint_float=float(solints_interval[s].replace('s','')) + for i in range(len(vislist)): + if solint_settings[solint]['interval'] != 'inf': + print(solint_settings[solint]['interval'][s]) + solint_float=float(solint_settings[solint]['interval'].replace('s','')) #use length of EB for S/N if inf_EB solint inf_EB_tint=solint_float - if solint_float > selfcal_library[selfcal_library['vislist'][i]]['TOS']: - inf_EB_tint=selfcal_library[selfcal_library['vislist'][i]]['TOS'] + if solint_float > selfcal_library[vislist[i]]['TOS']: + inf_EB_tint=selfcal_library[vislist[i]]['TOS'] else: - inf_EB_tint=selfcal_library[selfcal_library['vislist'][i]]['TOS'] + inf_EB_tint=selfcal_library[vislist[i]]['TOS'] SNR_self_EB[i]=SNR/((n_ant)**0.5*(selfcal_library['Total_TOS']/inf_EB_tint)**0.5) - SNR_self_EB_spw[selfcal_library['vislist'][i]]={} - SNR_self_EB_bb[selfcal_library['vislist'][i]]={} + SNR_self_EB_spw[vislist[i]]={} + SNR_self_EB_bb[vislist[i]]={} for spw in selfcal_library['spw_map']: - if selfcal_library['vislist'][i] in selfcal_library['spw_map'][spw]: - SNR_self_EB_spw[selfcal_library['vislist'][i]][str(spw)]=(polscale)**-0.5*SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/inf_EB_tint)**0.5)*(selfcal_library[selfcal_library['vislist'][i]]['per_spw_stats'][selfcal_library['spw_map'][spw][selfcal_library['vislist'][i]]]['effective_bandwidth']/selfcal_library[selfcal_library['vislist'][i]]['total_effective_bandwidth'])**0.5 - print(selfcal_library[selfcal_library['vislist'][i]]['baseband']) + if vislist[i] in selfcal_library['spw_map'][spw]: + SNR_self_EB_spw[vislist[i]][str(spw)]=(polscale)**-0.5*SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/inf_EB_tint)**0.5)*(selfcal_library[vislist[i]]['per_spw_stats'][selfcal_library['spw_map'][spw][vislist[i]]]['effective_bandwidth']/selfcal_library[vislist[i]]['total_effective_bandwidth'])**0.5 + print(selfcal_library[vislist[i]]['baseband']) print('SNR_self_EB_spw: ',SNR_self_EB_spw) - for baseband in selfcal_library[selfcal_library['vislist'][i]]['baseband']: - SNR_self_EB_bb[selfcal_library['vislist'][i]][baseband]=(polscale)**-0.5*SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/inf_EB_tint)**0.5)*(selfcal_library[selfcal_library['vislist'][i]]['baseband'][baseband]['total_effective_bandwidth']/selfcal_library[selfcal_library['vislist'][i]]['total_effective_bandwidth'])**0.5 + for baseband in selfcal_library[vislist[i]]['baseband']: + SNR_self_EB_bb[vislist[i]][baseband]=(polscale)**-0.5*SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/inf_EB_tint)**0.5)*(selfcal_library[vislist[i]]['baseband'][baseband]['total_effective_bandwidth']/selfcal_library[vislist[i]]['total_effective_bandwidth'])**0.5 print('SNR_self_EB_bb: ',SNR_self_EB_bb) for spw in selfcal_library['spw_map']: mean_SNR_spw=0.0 total_vis = 0 - for j in range(len(selfcal_library['vislist'])): - if selfcal_library['vislist'][j] in selfcal_library['spw_map'][spw]: - mean_SNR_spw+=SNR_self_EB_spw[selfcal_library['vislist'][j]][str(spw)] - total_vis += 1 - mean_SNR_spw=mean_SNR_spw/total_vis - solint_snr_per_spw[solint][str(spw)]=mean_SNR_spw + for j in range(len(vislist)): + if vislist[j] in selfcal_library['spw_map'][spw]: + mean_SNR_spw+=SNR_self_EB_spw[vislist[j]][str(spw)] + total_vis += 1 + if total_vis > 0: + mean_SNR_spw=mean_SNR_spw/total_vis + solint_snr_per_spw[solint][str(spw)]=mean_SNR_spw + else: + mean_SNR_spw = 0.0 for baseband in selfcal_library[selfcal_library['vislist'][i]]['baseband']: mean_SNR_bb=0.0 - for j in range(len(selfcal_library['vislist'])): - if baseband in SNR_self_EB_bb[selfcal_library['vislist'][j]].keys(): - mean_SNR_bb+=SNR_self_EB_bb[selfcal_library['vislist'][j]][baseband] - mean_SNR_bb=mean_SNR_bb/len(selfcal_library['vislist']) + for j in range(len(vislist)): + if baseband in SNR_self_EB_bb[vislist[j]].keys(): + mean_SNR_bb+=SNR_self_EB_bb[vislist[j]][baseband] + if len(vislist) > 0: + mean_SNR_bb=mean_SNR_bb/len(vislist) + else: + mean_SNR_bb = 0.0 print('mean_SNR_bb',mean_SNR_bb,baseband) solint_snr_per_bb[solint][baseband]=mean_SNR_bb solint_snr[solint]=np.mean(SNR_self_EB) selfcal_library['per_EB_SNR']=np.mean(SNR_self_EB) - elif solint =='scan_inf': - selfcal_library['per_scan_SNR']=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/selfcal_library['Median_scan_time'])**0.5) + elif solint_settings[solint]['sub-name'] =='scan_inf': + selfcal_library['per_scan_SNR']=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/selfcal_library[vislist[0]]['Median_scan_time'])**0.5) solint_snr[solint]=selfcal_library['per_scan_SNR'] for spw in selfcal_library['spw_map']: - vis = list(selfcal_library['spw_map'][spw].keys())[0] - true_spw = selfcal_library['spw_map'][spw][vis] - solint_snr_per_spw[solint][str(spw)]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/selfcal_library['Median_scan_time'])**0.5)*(selfcal_library[vis]['per_spw_stats'][true_spw]['effective_bandwidth']/selfcal_library[vis]['total_effective_bandwidth'])**0.5 + #vis = list(selfcal_library['spw_map'][spw].keys())[0] + vis = vislist[0] + if vis in selfcal_library['spw_map'][spw]: + true_spw = selfcal_library['spw_map'][spw][vis] + solint_snr_per_spw[solint][str(spw)]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/selfcal_library[vis]['Median_scan_time'])**0.5)*(selfcal_library[vis]['per_spw_stats'][true_spw]['effective_bandwidth']/selfcal_library[vis]['total_effective_bandwidth'])**0.5 for baseband in selfcal_library[vis]['baseband']: - solint_snr_per_bb[solint][baseband]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/selfcal_library['Median_scan_time'])**0.5)*(selfcal_library[vis]['baseband'][baseband]['total_effective_bandwidth']/selfcal_library[vis]['total_effective_bandwidth'])**0.5 - elif solint =='inf' or solint == 'inf_ap' or solint == 'inf_delay': - selfcal_library['per_scan_SNR']=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/(selfcal_library['Median_scan_time']/selfcal_library['Median_fields_per_scan']))**0.5) + solint_snr_per_bb[solint][baseband]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/selfcal_library[vis]['Median_scan_time'])**0.5)*(selfcal_library[vis]['baseband'][baseband]['total_effective_bandwidth']/selfcal_library[vis]['total_effective_bandwidth'])**0.5 + elif solint_settings[solint]['interval'] =='inf' or solint_settings[solint]['interval'] == 'inf_ap': + selfcal_library['per_scan_SNR']=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/(selfcal_library[vislist[0]]['Median_scan_time']/selfcal_library[vislist[0]]['Median_fields_per_scan']))**0.5) solint_snr[solint]=selfcal_library['per_scan_SNR'] for spw in selfcal_library['spw_map']: - vis = list(selfcal_library['spw_map'][spw].keys())[0] - true_spw = selfcal_library['spw_map'][spw][vis] - solint_snr_per_spw[solint][str(spw)]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/(selfcal_library['Median_scan_time']/selfcal_library['Median_fields_per_scan']))**0.5)*(selfcal_library[vis]['per_spw_stats'][true_spw]['effective_bandwidth']/selfcal_library[vis]['total_effective_bandwidth'])**0.5 + #vis = list(selfcal_library['spw_map'][spw].keys())[0] + vis = vislist[0] + if vis in selfcal_library['spw_map'][spw]: + true_spw = selfcal_library['spw_map'][spw][vis] + solint_snr_per_spw[solint][str(spw)]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/(selfcal_library[vis]['Median_scan_time']/selfcal_library[vis]['Median_fields_per_scan']))**0.5)*(selfcal_library[vis]['per_spw_stats'][true_spw]['effective_bandwidth']/selfcal_library[vis]['total_effective_bandwidth'])**0.5 for baseband in selfcal_library[vis]['baseband']: - solint_snr_per_bb[solint][baseband]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/(selfcal_library['Median_scan_time']/selfcal_library['Median_fields_per_scan']))**0.5)*(selfcal_library[vis]['baseband'][baseband]['total_effective_bandwidth']/selfcal_library[vis]['total_effective_bandwidth'])**0.5 - elif solint == 'int' or solint == 'int_ap': + solint_snr_per_bb[solint][baseband]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/(selfcal_library[vis]['Median_scan_time']/selfcal_library[vis]['Median_fields_per_scan']))**0.5)*(selfcal_library[vis]['baseband'][baseband]['total_effective_bandwidth']/selfcal_library[vis]['total_effective_bandwidth'])**0.5 + elif solint_settings[solint]['interval'] == 'int': solint_snr[solint]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/integration_time)**0.5) for spw in selfcal_library['spw_map']: - vis = list(selfcal_library['spw_map'][spw].keys())[0] - true_spw = selfcal_library['spw_map'][spw][vis] - solint_snr_per_spw[solint][str(spw)]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/integration_time)**0.5)*(selfcal_library[vis]['per_spw_stats'][true_spw]['effective_bandwidth']/selfcal_library[vis]['total_effective_bandwidth'])**0.5 + #vis = list(selfcal_library['spw_map'][spw].keys())[0] + vis = vislist[0] + if vis in selfcal_library['spw_map'][spw]: + true_spw = selfcal_library['spw_map'][spw][vis] + solint_snr_per_spw[solint][str(spw)]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/integration_time)**0.5)*(selfcal_library[vis]['per_spw_stats'][true_spw]['effective_bandwidth']/selfcal_library[vis]['total_effective_bandwidth'])**0.5 for baseband in selfcal_library[vis]['baseband']: solint_snr_per_bb[solint][baseband]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/integration_time)**0.5)*(selfcal_library[vis]['baseband'][baseband]['total_effective_bandwidth']/selfcal_library[vis]['total_effective_bandwidth'])**0.5 else: - solint_float=float(solint.replace('s','').replace('_ap','')) + solint_float=float(solint_settings[solint]['interval'].replace('s','').replace('_ap','')) solint_snr[solint]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/solint_float)**0.5) for spw in selfcal_library['spw_map']: - vis = list(selfcal_library['spw_map'][spw].keys())[0] - true_spw = selfcal_library['spw_map'][spw][vis] - solint_snr_per_spw[solint][str(spw)]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/solint_float)**0.5)*(selfcal_library[vis]['per_spw_stats'][true_spw]['effective_bandwidth']/selfcal_library[vis]['total_effective_bandwidth'])**0.5 + #vis = list(selfcal_library['spw_map'][spw].keys())[0] + vis = vislist[0] + if vis in selfcal_library['spw_map'][spw]: + true_spw = selfcal_library['spw_map'][spw][vis] + solint_snr_per_spw[solint][str(spw)]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/solint_float)**0.5)*(selfcal_library[vis]['per_spw_stats'][true_spw]['effective_bandwidth']/selfcal_library[vis]['total_effective_bandwidth'])**0.5 for baseband in selfcal_library[vis]['baseband']: solint_snr_per_bb[solint][baseband]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/solint_float)**0.5)*(selfcal_library[vis]['baseband'][baseband]['total_effective_bandwidth']/selfcal_library[vis]['total_effective_bandwidth'])**0.5 return solint_snr,solint_snr_per_spw,solint_snr_per_bb -def get_SNR_self_update(selfcal_library,selfcal_plan,n_ant,solint_curr,solint_next,integration_time,solint_snr): +def get_SNR_self_update(vis,selfcal_library,selfcal_plan,n_ant,solint_curr,solint_next,integration_time,solint_snr): + """ maxspws=0 maxspwvis='' for vis in selfcal_library['vislist']: if selfcal_library[vis]['n_spws'] >= maxspws: maxspws=selfcal_library[vis]['n_spws'] maxspwvis=vis+'' - SNR = max(selfcal_library[selfcal_library['vislist'][0]][solint_curr]['SNR_post'],selfcal_library[selfcal_library['vislist'][0]][solint_curr]['intflux_post']/selfcal_library[selfcal_library['vislist'][0]][solint_curr]['e_intflux_post']) + """ + SNR = max(selfcal_library[vis][solint_curr]['SNR_post'],selfcal_library[vis][solint_curr]['intflux_post']/selfcal_library[vis][solint_curr]['e_intflux_post']) SNR_orig = max(selfcal_library['SNR_orig'],selfcal_library['intflux_orig']/selfcal_library['e_intflux_orig']) SNR_ratio = SNR / SNR_orig #solint_snr[solint_next]=SNR_ratio*solint_snr[solint_next] + #for vis in selfcal_library['vislist']: solint_snr[solint_next]=SNR_ratio*solint_snr[solint_next] for spw in selfcal_library['spw_map']: - selfcal_plan['solint_snr_per_spw'][solint_next][str(spw)]=selfcal_plan['solint_snr_per_spw'][solint_next][str(spw)]*SNR_ratio + if vis in selfcal_library['spw_map'][spw]: + selfcal_plan[vis]['solint_snr_per_spw'][solint_next][str(spw)]=selfcal_plan[vis]['solint_snr_per_spw'][solint_next][str(spw)]*SNR_ratio for baseband in selfcal_library[vis]['baseband']: - selfcal_plan['solint_snr_per_bb'][solint_next][baseband]=selfcal_plan['solint_snr_per_bb'][solint_next][baseband]*SNR_ratio + selfcal_plan[vis]['solint_snr_per_bb'][solint_next][baseband]=selfcal_plan[vis]['solint_snr_per_bb'][solint_next][baseband]*SNR_ratio def get_sensitivity(vislist,selfcal_library,field='',virtual_spw='all',chan=0,cellsize='0.025arcsec',imsize=[1600,1600],robust=0.5,specmode='mfs',uvtaper=''): @@ -3771,7 +3797,7 @@ def select_best_gaincal_mode(selfcal_library,selfcal_plan,vis,gaintable_prefix,s # Check whether any spws have estimated SNR < 3, in which case we should not (initially) allow 'per_spw' coarsest_solint=selfcal_plan['solints'][0] # use this instead of assuming inf_EB - if preferred_mode == 'per_spw' and np.any([selfcal_plan['solint_snr_per_spw'][coarsest_solint][str(selfcal_library['reverse_spw_map'][vis][int(spw)])] < \ + if preferred_mode == 'per_spw' and np.any([selfcal_plan[vis]['solint_snr_per_spw'][coarsest_solint][str(selfcal_library['reverse_spw_map'][vis][int(spw)])] < \ minsnr_to_proceed for spw in spwlist]): if 'per_bb' in selfcal_plan[vis]['solint_settings'][solint]['modes_to_attempt']: preferred_mode = 'per_bb' @@ -3786,7 +3812,7 @@ def select_best_gaincal_mode(selfcal_library,selfcal_plan,vis,gaintable_prefix,s for i in range(len(spwlist)): # use >= to not always map if an spw has flagged solutions for a given antenna if np.min(selfcal_plan[vis]['solint_settings'][solint]['delta_nflags']['per_spw'][i]) >= max_flagged_ants_spwmap or \ - selfcal_plan['solint_snr_per_spw'][coarsest_solint][str(selfcal_library['reverse_spw_map'][vis][int(spwlist[i])])] < minsnr_to_proceed or \ + selfcal_plan[vis]['solint_snr_per_spw'][coarsest_solint][str(selfcal_library['reverse_spw_map'][vis][int(spwlist[i])])] < minsnr_to_proceed or \ selfcal_plan[vis]['solint_settings'][solint]['fracflagged']['per_spw'][i] == 1.0: fallback='spwmap' spwmap[i]=1.0 @@ -5004,9 +5030,11 @@ def get_min_SNR_spw(snr_per_spw): def remove_modes(selfcal_plan,vis,start_index): # remove the per_spw and/or per_bb modes for solints following current solint preferred_mode=selfcal_plan[vis]['solint_settings'][selfcal_plan['solints'][start_index]]['final_mode'] for j in range(start_index+1,len(selfcal_plan['solints'])): + if selfcal_plan['solints'][j] not in selfcal_plan[vis]['solint_settings']: + continue if 'ap' in selfcal_plan['solints'][j] and 'ap' not in selfcal_plan['solints'][start_index]: # exempt over ap solints since they go back to a longer solint continue - if 'delay' in selfcal_plan['solints'][j]: # do not remove for delay solints since they cannot use combinespw + if 'd' in selfcal_plan['solints'][j]: # do not remove for delay solints since they cannot use combinespw continue if preferred_mode == 'per_bb' or preferred_mode == 'combinespw': if 'per_spw' in selfcal_plan[vis]['solint_settings'][selfcal_plan['solints'][j]]['modes_to_attempt']: diff --git a/auto_selfcal/tests/test_auto_selfcal.py b/auto_selfcal/tests/test_auto_selfcal.py index 8910247..f4737dc 100644 --- a/auto_selfcal/tests/test_auto_selfcal.py +++ b/auto_selfcal/tests/test_auto_selfcal.py @@ -36,7 +36,6 @@ pytest.param('disjoint_missing_spwfid_combo', id='disjoint_missing_spwfid_combo'), # tests spw-mapping, particularly in cases where spws and sub-fields are missing pytest.param('EE11.1.00032.S', id='EE11.1.00032.S'), # tests a dataset with spw names crossing a naming scheme change pytest.param('VLBA_CASA_Guide', id='VLBA_CASA_Guide'), # VLBA dataset - ] ) def test_benchmark(tmp_path, dataset): @@ -62,7 +61,7 @@ def test_benchmark(tmp_path, dataset): if 'VLBA' in dataset: job = ex.submit(auto_selfcal, sort_targets_and_EBs=True, weblog=True, parallel=True, iscalibrator=True, targets='J1154+6022',applytargets='J1203+6031',imsize=640,cell='0.0002arcsec') else: - job = ex.submit(auto_selfcal, sort_targets_and_EBs=True, weblog=True, parallel=True) + job = ex.submit(auto_selfcal, sort_targets_and_EBs=True, weblog=True, parallel=True, uniform_solints=True) job.wait() assert job.state in ['DONE','COMPLETED'] @@ -72,10 +71,28 @@ def test_benchmark(tmp_path, dataset): with open('selfcal_library.pickle', 'rb') as handle: selfcal_library2 = pickle.load(handle) - difference_count = compare_two_dictionaries(selfcal_library1, selfcal_library2, tolerance=1e-3,\ + with open('selfcal_plan.pickle', 'rb') as handle: + selfcal_plan = pickle.load(handle) + + solint_map = {} + for target in selfcal_library2: + for band in selfcal_library2[target]: + for vis in selfcal_library2[target][band]['vislist']: + for solint in selfcal_plan[target][band][vis]['solint_settings']: + if solint not in solint_map: + solint_map[solint] = [] + + mapped_solint = selfcal_plan[target][band][vis]['solint_settings'][solint]['sub-name'] + + solint_map[solint].append(mapped_solint) + print(solint_map) + + difference_count = compare_two_dictionaries(selfcal_library1, selfcal_library2, tolerance=1e-3, key_map=solint_map,\ exclude=['vislist_orig','field_str','imsize','flux_threshold','overlap_tol','bands_for_targets',\ 'am_dogrowprune','am_growiterations','am_lownoisethreshold','am_minbeamfrac',\ - 'am_noisethreshold','am_sidelobethreshold','am_smoothfactor','telescope']) + 'am_noisethreshold','am_sidelobethreshold','am_smoothfactor','telescope', "final_phase_solint", "final_solint", \ + "gaintable_final", "per_EB_SNR", "per_scan_SNR", "vislist-to-gaincal", "telescope", "gaintable","sub-fields-to-gaincal", \ + "sub-fields-to-selfcal", "Stop_Reason"]) for msfile in starting_MS_files: print(f"rm -rf {os.path.basename(msfile)}") @@ -131,10 +148,27 @@ def test_on_github(tmp_path, request, zip_file, link): with open('selfcal_library.pickle', 'rb') as handle: selfcal_library2 = pickle.load(handle) - difference_count = compare_two_dictionaries(selfcal_library1, selfcal_library2, tolerance=0.001,\ - exclude=['vislist_orig','field_str','imsize','flux_threshold','overlap_tol','bands_for_targets',\ - 'am_dogrowprune','am_growiterations','am_lownoisethreshold','am_minbeamfrac',\ - 'am_noisethreshold','am_sidelobethreshold','am_smoothfactor','telescope']) + with open('selfcal_plan.pickle', 'rb') as handle: + selfcal_plan = pickle.load(handle) + + solint_map = {} + for target in selfcal_library2: + for band in selfcal_library2[target]: + for vis in selfcal_library2[target][band]['vislist']: + for solint in selfcal_plan[target][band][vis]['solint_settings']: + if solint not in solint_map: + solint_map[solint] = [] + + mapped_solint = selfcal_plan[target][band][vis]['solint_settings'][solint]['sub-name'] + + solint_map[solint].append(mapped_solint) + #print(solint_map) + + difference_count = compare_two_dictionaries(selfcal_library1, selfcal_library2, tolerance=0.001, key_map=solint_map, + exclude=["final_phase_solint", "final_solint", "gaintable_final", "per_EB_SNR", "vislist-to-gaincal", "telescope", + "gaintable","sub-fields-to-gaincal", "sub-fields-to-selfcal", "am_dogrowprune", "am_growiterations", + "am_lownoisethreshold", "am_minbeamfrac", "am_noisethreshold", "am_sidelobethreshold", + "am_smoothfactor", "Stop_Reason"]) assert difference_count == 0 @@ -150,13 +184,15 @@ def compare_values(list1, list2, tol=1e-3): return np.all([compare_values(list1[i], list2[i], tol=tol) for i in range(len(list1))]) elif type(list1) == str or type(list1) == np.str_ or type(list1) == bool: return list1 == list2 + elif type(list1) == type(None): + return list1 == list2 else: if list1 == 0: return abs(list2) < tol else: return abs(list1 - list2) < abs(list1*tol) -def compare_two_dictionaries(dictionary1, dictionary2, path=[], exclude=[], tolerance=1e-3): +def compare_two_dictionaries(dictionary1, dictionary2, path=[], exclude=[], tolerance=1e-3, key_map={}, verbose=True): if isinstance(dictionary1, str): with open(dictionary1, 'rb') as handle: dictionary1 = pickle.load(handle) @@ -172,35 +208,64 @@ def compare_two_dictionaries(dictionary1, dictionary2, path=[], exclude=[], tole if key in exclude: continue - if key not in intersect_keys: + if key not in intersect_keys and key not in key_map and not np.any([key in key_map[k] for k in key_map]): if key not in dictionary1: - print('/'.join([str(p) for p in path])+"/"+key+" not in dictionary1") + if verbose: + print('/'.join([str(p) for p in path])+"/"+key+" not in dictionary1") else: - print('/'.join([str(p) for p in path])+"/"+key+" not in dictionary2") + if verbose: + print('/'.join([str(p) for p in path])+"/"+key+" not in dictionary2") difference_count += 1 continue + elif key not in intersect_keys and key not in key_map and np.any([key in key_map[k] for k in key_map]): + #print(f'key {key} has changed in dictionary2 and will be matched elsewhere') + continue - if key not in dictionary1 and int(key) in dictionary1: - key = int(key) + try: + if key not in dictionary1 and key not in key_map and int(key) in dictionary1: + key = int(key) + except: + continue - if type(dictionary1[key]) == dict: - difference_count += compare_two_dictionaries(dictionary1[key], dictionary2[key], path.copy()+[key], exclude=exclude, tolerance=tolerance) + if key in dictionary2 and not key in dictionary1 and key in key_map: + #print(f'Checking whether key {key} has its name changed') + found = False + for alt_key in key_map[key]: + #print(f'Checking for {alt_key} in dictionary1') + if alt_key in dictionary1: + found = True + break + + if found: + pass + #print(f"Using alternative key {alt_key} to match with key {key}") + else: + if verbose: + print(f"No match for key {key} found in dictionary1, this is a difference") + difference_count += 1 + continue + else: + alt_key = key + + if type(dictionary1[alt_key]) == dict: + difference_count += compare_two_dictionaries(dictionary1[alt_key], dictionary2[key], path.copy()+[key], exclude=exclude, tolerance=tolerance, key_map=key_map, verbose=verbose) else: - value1 = np.array(dictionary1[key])[np.argsort(dictionary1['vislist'])] if key in ['spws_per_vis','vislist'] else dictionary1[key] + value1 = np.array(dictionary1[alt_key])[np.argsort(dictionary1['vislist'])] if alt_key in ['spws_per_vis','vislist'] else dictionary1[alt_key] value2 = np.array(dictionary2[key])[np.argsort(dictionary2['vislist'])] if key in ['spws_per_vis','vislist'] else dictionary2[key] #value1 = np.array(dictionary1[key])[np.argsort(dictionary1['vislist'])] if key in ['spws_per_vis'] else dictionary1[key] #value2 = np.array(dictionary2[key])[np.argsort(dictionary2['vislist'])] if key in ['spws_per_vis'] else dictionary2[key] if key == 'gaincal_combine': - value1 = dictionary1[key].split(',') + value1 = dictionary1[alt_key].split(',') value1.sort() value2 = dictionary2[key].split(',') value2.sort() if not compare_values(value1, value2, tol=tolerance): - print('/'.join([str(p) for p in path])+"/"+key, dictionary1[key], dictionary2[key]) + if verbose: + print('/'.join([str(p) for p in path])+"/"+key, dictionary1[key], dictionary2[key]) difference_count += 1 return difference_count diff --git a/auto_selfcal/weblog_creation.py b/auto_selfcal/weblog_creation.py index 19e3fff..7debcbb 100644 --- a/auto_selfcal/weblog_creation.py +++ b/auto_selfcal/weblog_creation.py @@ -14,7 +14,7 @@ from casatools import msmetadata as msmdtool from casatools import table as tbtool from casatools import ms as mstool -from casaviewer import imview +#from casaviewer import imview from PIL import Image ms = mstool() @@ -79,15 +79,17 @@ def generate_weblog(sclib,selfcal_plan,directory='weblog'): htmlOut.writelines('Aligned EBs?: False\n') htmlOut.writelines('Selfcal Success?: '+str(sclib[target][band]['SC_success'])+'
\n') keylist=sclib[target][band].keys() - if 'Stop_Reason' not in keylist: + if np.all(['Stop_Reason' not in sclib[target][band][vis].keys() for vis in sclib[target][band]['vislist']]): htmlOut.writelines('Stop Reason: Estimated Selfcal S/N too low for solint

\n') if sclib[target][band]['SC_success']==False: render_summary_table(htmlOut,sclib,target,band,directory=directory) continue else: - htmlOut.writelines('Stop Reason: '+str(sclib[target][band]['Stop_Reason'])+'

\n') - print(target,band,sclib[target][band]['Stop_Reason']) - if (('Estimated_SNR_too_low_for_solint' in sclib[target][band]['Stop_Reason']) or ('Selfcal_Not_Attempted' in sclib[target][band]['Stop_Reason'])) and sclib[target][band]['final_solint']=='None': + for vis in sclib[target][band]['vislist']: + htmlOut.writelines(vis + 'Stop Reason: '+str(sclib[target][band][vis]['Stop_Reason'])+'

\n') + print(vis, target,band,sclib[target][band][vis]['Stop_Reason']) + if ((np.all(['Estimated_SNR_too_low_for_solint' in sclib[target][band][vis]['Stop_Reason'] for vis in sclib[target][band]['vislist']])) or \ + (np.all(['Selfcal_Not_Attempted' in sclib[target][band][vis]['Stop_Reason'] for vis in sclib[target][band]['vislist']]))) and sclib[target][band]['final_solint']=='None': render_summary_table(htmlOut,sclib,target,band,directory=directory) continue htmlOut.writelines('Final Successful solint: '+str(sclib[target][band]['final_solint'])+'

\n') @@ -115,7 +117,7 @@ def generate_weblog(sclib,selfcal_plan,directory='weblog'): htmlOut.writelines('Noise Characteristics
\n') # Solint summary table - if 'Empty model' not in sclib[target][band]['Stop_Reason']: + if 'Empty model' not in sclib[target][band][sclib[target][band]['vislist'][0]]['Stop_Reason']: render_selfcal_solint_summary_table(htmlOut,sclib,target,band,selfcal_plan) # PER SPW STATS TABLE @@ -128,7 +130,7 @@ def generate_weblog(sclib,selfcal_plan,directory='weblog'): htmlOut.close() # Pages for each solint - if 'Empty model' not in sclib[target][band]['Stop_Reason']: + if 'Empty model' not in sclib[target][band][sclib[target][band]['vislist'][0]]['Stop_Reason']: render_per_solint_QA_pages(sclib,selfcal_plan,bands,directory=directory) @@ -238,8 +240,36 @@ def render_selfcal_solint_summary_table(htmlOut,sclib,target,band,selfcal_plan): line+=''+solint+'\n ' line+='\n' htmlOut.writelines(line) - vis_keys=list(sclib[target][band][vislist[len(vislist)-1]].keys()) - quantities=['Pass','intflux_final','intflux_improvement','SNR_final','SNR_Improvement','SNR_NF_final','SNR_NF_Improvement','RMS_final','RMS_Improvement','RMS_NF_final','RMS_NF_Improvement','Beam_Ratio','clean_threshold','Plots'] + htmlOut.writelines('\n Solution interval by EB: \n') + for vis in vislist: + line=f'\n {vis}: \n' + for solint in solint_list: + if solint in selfcal_plan[target][band][vis]['solint_settings']: + line += f' {selfcal_plan[target][band][vis]["solint_settings"][solint]["interval"]} \n' + else: + line += ' - \n' + line += '\n' + htmlOut.writelines(line) + htmlOut.writelines('\n Result: \n') + for vis in vislist: + line=f'\n {vis}: \n' + for solint in solint_list: + if solint in sclib[target][band][vis]: + if 'Pass' in sclib[target][band][vis][solint]: + if sclib[target][band][vis][solint]['Pass'] == False: + line+=' {} {}\n'.format('Fail',sclib[target][band][vis][solint]['Fail_Reason']) + elif sclib[target][band][vis][solint]['Pass'] == 'None': + line+=' {} {}\n'.format('Not attempted',sclib[target][band][vis][solint]['Fail_Reason']) + else: + line+=' {}\n'.format('Pass') + else: + line+=' {}\n'.format('None') + else: + line += ' - \n' + line += '\n' + htmlOut.writelines(line) + htmlOut.writelines('\n Selfcal stats: \n') + quantities=['intflux_final','intflux_improvement','SNR_final','SNR_Improvement','SNR_NF_final','SNR_NF_Improvement','RMS_final','RMS_Improvement','RMS_NF_final','RMS_NF_Improvement','Beam_Ratio','clean_threshold','Plots'] for key in quantities: if key =='Pass': line='\n Result: \n' @@ -270,51 +300,61 @@ def render_selfcal_solint_summary_table(htmlOut,sclib,target,band,selfcal_plan): if key =='Plots': line='\n Plots: \n' for solint in solint_list: + if np.any([solint in sclib[target][band][vis] for vis in vislist]): + if np.any([sclib[target][band][vis][solint]['Pass'] != 'None' for vis in vislist]): + ivis = np.where([solint in sclib[target][band][vis] and sclib[target][band][vis][solint]['Pass'] != 'None' for vis in vislist])[0][0] + else: + ivis = np.where([solint in sclib[target][band][vis] for vis in vislist])[0][0] + else: + ivis = len(vislist)-1 + + vis_keys=list(sclib[target][band][vislist[ivis]].keys()) + if solint in vis_keys: - vis_solint_keys=sclib[target][band][vislist[len(vislist)-1]][solint].keys() - if key != 'Pass' and sclib[target][band][vislist[len(vislist)-1]][solint]['Pass'] == 'None': + vis_solint_keys=sclib[target][band][vislist[ivis]][solint].keys() + if key != 'Pass' and sclib[target][band][vislist[ivis]][solint]['Pass'] == 'None': line+=' - \n' continue if key=='Pass': - if key in sclib[target][band][vislist[len(vislist)-1]][solint]: - if sclib[target][band][vislist[len(vislist)-1]][solint]['Pass'] == False: - line+=' {} {}\n'.format('Fail',sclib[target][band][vislist[len(vislist)-1]][solint]['Fail_Reason']) - elif sclib[target][band][vislist[len(vislist)-1]][solint]['Pass'] == 'None': - line+=' {} {}\n'.format('Not attempted',sclib[target][band][vislist[len(vislist)-1]][solint]['Fail_Reason']) + if key in sclib[target][band][vislist[ivis]][solint]: + if sclib[target][band][vislist[ivis]][solint]['Pass'] == False: + line+=' {} {}\n'.format('Fail',sclib[target][band][vislist[ivis]][solint]['Fail_Reason']) + elif sclib[target][band][vislist[ivis]][solint]['Pass'] == 'None': + line+=' {} {}\n'.format('Not attempted',sclib[target][band][vislist[ivis]][solint]['Fail_Reason']) else: line+=' {}\n'.format('Pass') else: line+=' {}\n'.format('None') if key=='intflux_final': - line+=' {:0.3f} +/- {:0.3f} mJy\n'.format(sclib[target][band][vislist[len(vislist)-1]][solint]['intflux_post']*1000.0,sclib[target][band][vislist[len(vislist)-1]][solint]['e_intflux_post']*1000.0) + line+=' {:0.3f} +/- {:0.3f} mJy\n'.format(sclib[target][band][vislist[ivis]][solint]['intflux_post']*1000.0,sclib[target][band][vislist[ivis]][solint]['e_intflux_post']*1000.0) if key=='intflux_improvement': - if sclib[target][band][vislist[len(vislist)-1]][solint]['intflux_pre'] == 0: + if sclib[target][band][vislist[ivis]][solint]['intflux_pre'] == 0: line+=' {:0.3f}\n'.format(1.0) else: - line+=' {:0.3f}\n'.format(sclib[target][band][vislist[len(vislist)-1]][solint]['intflux_post']/sclib[target][band][vislist[len(vislist)-1]][solint]['intflux_pre']) + line+=' {:0.3f}\n'.format(sclib[target][band][vislist[ivis]][solint]['intflux_post']/sclib[target][band][vislist[ivis]][solint]['intflux_pre']) if key=='SNR_final': - line+=' {:0.3f}\n'.format(sclib[target][band][vislist[len(vislist)-1]][solint]['SNR_post']) + line+=' {:0.3f}\n'.format(sclib[target][band][vislist[ivis]][solint]['SNR_post']) if key=='SNR_Improvement': - line+=' {:0.3f}\n'.format(sclib[target][band][vislist[len(vislist)-1]][solint]['SNR_post']/sclib[target][band][vislist[len(vislist)-1]][solint]['SNR_pre']) + line+=' {:0.3f}\n'.format(sclib[target][band][vislist[ivis]][solint]['SNR_post']/sclib[target][band][vislist[ivis]][solint]['SNR_pre']) if key=='SNR_NF_final': - line+=' {:0.3f}\n'.format(sclib[target][band][vislist[len(vislist)-1]][solint]['SNR_NF_post']) + line+=' {:0.3f}\n'.format(sclib[target][band][vislist[ivis]][solint]['SNR_NF_post']) if key=='SNR_NF_Improvement': - line+=' {:0.3f}\n'.format(sclib[target][band][vislist[len(vislist)-1]][solint]['SNR_NF_post']/sclib[target][band][vislist[len(vislist)-1]][solint]['SNR_NF_pre']) + line+=' {:0.3f}\n'.format(sclib[target][band][vislist[ivis]][solint]['SNR_NF_post']/sclib[target][band][vislist[ivis]][solint]['SNR_NF_pre']) if key=='RMS_final': - line+=' {:0.3e} mJy/bm\n'.format(sclib[target][band][vislist[len(vislist)-1]][solint]['RMS_post']*1000.0) + line+=' {:0.3e} mJy/bm\n'.format(sclib[target][band][vislist[ivis]][solint]['RMS_post']*1000.0) if key=='RMS_Improvement': - line+=' {:0.3e}\n'.format(sclib[target][band][vislist[len(vislist)-1]][solint]['RMS_pre']/sclib[target][band][vislist[len(vislist)-1]][solint]['RMS_post']) + line+=' {:0.3e}\n'.format(sclib[target][band][vislist[ivis]][solint]['RMS_pre']/sclib[target][band][vislist[ivis]][solint]['RMS_post']) if key=='RMS_NF_final': - line+=' {:0.3e} mJy/bm\n'.format(sclib[target][band][vislist[len(vislist)-1]][solint]['RMS_NF_post']*1000.0) + line+=' {:0.3e} mJy/bm\n'.format(sclib[target][band][vislist[ivis]][solint]['RMS_NF_post']*1000.0) if key=='RMS_NF_Improvement': - line+=' {:0.3e}\n'.format(sclib[target][band][vislist[len(vislist)-1]][solint]['RMS_NF_pre']/sclib[target][band][vislist[len(vislist)-1]][solint]['RMS_NF_post']) + line+=' {:0.3e}\n'.format(sclib[target][band][vislist[ivis]][solint]['RMS_NF_pre']/sclib[target][band][vislist[ivis]][solint]['RMS_NF_post']) if key=='Beam_Ratio': - line+=' {:0.3e}\n'.format((sclib[target][band][vislist[len(vislist)-1]][solint]['Beam_major_post']*sclib[target][band][vislist[len(vislist)-1]][solint]['Beam_minor_post'])/(sclib[target][band]['Beam_major_orig']*sclib[target][band]['Beam_minor_orig'])) + line+=' {:0.3e}\n'.format((sclib[target][band][vislist[ivis]][solint]['Beam_major_post']*sclib[target][band][vislist[ivis]][solint]['Beam_minor_post'])/(sclib[target][band]['Beam_major_orig']*sclib[target][band]['Beam_minor_orig'])) if key =='clean_threshold': if key in vis_solint_keys: - line+=' {:0.3e} mJy/bm\n'.format(sclib[target][band][vislist[len(vislist)-1]][solint]['clean_threshold']*1000.0) + line+=' {:0.3e} mJy/bm\n'.format(sclib[target][band][vislist[ivis]][solint]['clean_threshold']*1000.0) else: line+=' Not Available\n' if key =='Plots': @@ -328,7 +368,7 @@ def render_selfcal_solint_summary_table(htmlOut,sclib,target,band,selfcal_plan): for vis in vislist: line='\n '+vis+': \n' for solint in solint_list: - if solint in vis_keys and sclib[target][band][vis][solint]['Pass'] != 'None' and 'gaintable' in sclib[target][band][vis][solint]: + if solint in sclib[target][band][vis] and sclib[target][band][vis][solint]['Pass'] != 'None' and 'gaintable' in sclib[target][band][vis][solint]: # only evaluate last gaintable not the pre-apply table gaintable=sclib[target][band][vis][solint]['gaintable'][len(sclib[target][band][vis][solint]['gaintable'])-1] line+='antenna positions with flagging plot\n' @@ -339,7 +379,7 @@ def render_selfcal_solint_summary_table(htmlOut,sclib,target,band,selfcal_plan): for quantity in ['Nsols_with_preflagged_data','Flagged_Sols_with_preflagged_data','Frac_Flagged_with_preflagged_data','Nsols_without_preflagged_data','Flagged_Sols_without_preflagged_data','Frac_Flagged_without_preflagged_data','Frac_sols_not_equiv_zero_per_spw','Frac_sols_not_equiv_zero_per_bb','SPW_Combine_Mode']: line='\n '+quantity+'\n' for solint in solint_list: - if solint in vis_keys and sclib[target][band][vis][solint]['Pass'] != 'None' and 'gaintable' in sclib[target][band][vis][solint]: + if solint in sclib[target][band][vis] and sclib[target][band][vis][solint]['Pass'] != 'None' and 'gaintable' in sclib[target][band][vis][solint]: # only evaluate last gaintable not the pre-apply table #gaintable=sclib[target][band][vis][solint]['gaintable'][len(sclib[target][band][vis][solint]['gaintable'])-1] #nflagged_sols, nsols=get_sols_flagged_solns(gaintable) @@ -483,7 +523,8 @@ def render_per_solint_QA_pages(sclib,selfcal_plan,bands,directory='weblog'): final_solint_to_plot=selfcal_plan[target][band]['solints'][final_solint_index+index_addition-1] - keylist=sclib[target][band][vislist[0]].keys() + #keylist=sclib[target][band][vislist[0]].keys() + keylist = [solint for solint in selfcal_plan[target][band]['solints'] if np.any([solint in sclib[target][band][vis] for vis in vislist])] if index_addition == 2 and final_solint_to_plot not in keylist: index_addition=index_addition-1 @@ -492,7 +533,10 @@ def render_per_solint_QA_pages(sclib,selfcal_plan,bands,directory='weblog'): #for i in range(final_solint_index+index_addition): for i in range(len(selfcal_plan[target][band]['solints'])): - if selfcal_plan[target][band]['solints'][i] not in keylist or sclib[target][band][vislist[len(vislist)-1]][selfcal_plan[target][band]['solints'][i]]['Pass'] == 'None': + representative_vislist = [vis for vis in vislist if selfcal_plan[target][band]['solints'][i] in sclib[target][band][vis]] + if len(representative_vislist) == 0: + continue + if selfcal_plan[target][band]['solints'][i] not in keylist or sclib[target][band][representative_vislist[0]][selfcal_plan[target][band]['solints'][i]]['Pass'] == 'None': continue htmlOutSolint=open(directory+'/'+target+'_'+band+'_'+selfcal_plan[target][band]['solints'][i]+'.html','w') htmlOutSolint.writelines('\n') @@ -504,22 +548,29 @@ def render_per_solint_QA_pages(sclib,selfcal_plan,bands,directory='weblog'): htmlOutSolint.writelines('

'+target+' Plots

\n') htmlOutSolint.writelines('

'+band+'

\n') htmlOutSolint.writelines('

Targets:

\n') - keylist=sclib[target][band][vislist[0]].keys() + #keylist=sclib[target][band][vislist[0]].keys() solints_string='' + print(keylist) for j in range(final_solint_index+index_addition): if selfcal_plan[target][band]['solints'][j] not in keylist: continue solints_string+=''+selfcal_plan[target][band]['solints'][j]+'
\n' htmlOutSolint.writelines('
Solints: '+solints_string) - htmlOutSolint.writelines('

Solint: '+selfcal_plan[target][band]['solints'][i]+'

\n') + htmlOutSolint.writelines('

Solint: '+selfcal_plan[target][band]['solints'][i]+'

\n') + for ivis, vis in enumerate(representative_vislist): + if ivis == len(representative_vislist)-1: + margin_string = 'style="margin-top: 0; padding-top:0;"' + else: + margin_string = 'style="margin : 0; padding-top:0;"' + htmlOutSolint.writelines(f'

{vis}: {selfcal_plan[target][band][vis]["solint_settings"][selfcal_plan[target][band]["solints"][i]]["interval"]}

\n') keylist_top=sclib[target][band].keys() htmlOutSolint.writelines('Back to Main Target/Band
\n') #must select last key for pre Jan 14th runs since they only wrote pass to the last MS dictionary entry - if "Pass" in sclib[target][band][vislist[len(vislist)-1]][selfcal_plan[target][band]['solints'][i]]: - passed=sclib[target][band][vislist[len(vislist)-1]][selfcal_plan[target][band]['solints'][i]]['Pass'] + if "Pass" in sclib[target][band][representative_vislist[-1]][selfcal_plan[target][band]['solints'][i]]: + passed=sclib[target][band][representative_vislist[-1]][selfcal_plan[target][band]['solints'][i]]['Pass'] else: passed = 'None' @@ -538,7 +589,7 @@ def render_per_solint_QA_pages(sclib,selfcal_plan,bands,directory='weblog'): htmlOutSolint.writelines('

Passed: True

\n') else: htmlOutSolint.writelines('

Passed: False

\n') - if 'Empty model' in sclib[target][band]['Stop_Reason']: + if 'Empty model' in sclib[target][band][sclib[target][band]['vislist'][0]]['Stop_Reason']: htmlOutSolint.writelines('Empty model image, no gains solved
\n') htmlOutSolint.writelines('\n') htmlOutSolint.writelines('\n') @@ -555,19 +606,19 @@ def render_per_solint_QA_pages(sclib,selfcal_plan,bands,directory='weblog'): htmlOutSolint.writelines('pre-SC-solint image\n') htmlOutSolint.writelines('pre-SC-solint image
\n') - htmlOutSolint.writelines('Post SC SNR: {:0.3f}'.format(sclib[target][band][vislist[0]][selfcal_plan[target][band]['solints'][i]]['SNR_post'])+'
Pre SC SNR: {:0.3f}'.format(sclib[target][band][vislist[0]][selfcal_plan[target][band]['solints'][i]]['SNR_pre'])+'

\n') - htmlOutSolint.writelines('Post SC RMS: {:0.7f}'.format(sclib[target][band][vislist[0]][selfcal_plan[target][band]['solints'][i]]['RMS_post'])+' Jy/beam
Pre SC RMS: {:0.7f}'.format(sclib[target][band][vislist[0]][selfcal_plan[target][band]['solints'][i]]['RMS_pre'])+' Jy/beam
\n') - htmlOutSolint.writelines('Post Beam: {:0.4f}"x{:0.4f}" {:0.3f} deg'.format(sclib[target][band][vislist[0]][selfcal_plan[target][band]['solints'][i]]['Beam_major_post'],sclib[target][band][vislist[0]][selfcal_plan[target][band]['solints'][i]]['Beam_minor_post'],sclib[target][band][vislist[0]][selfcal_plan[target][band]['solints'][i]]['Beam_PA_post'])+'
\n') - htmlOutSolint.writelines('Pre Beam: {:0.4f}"x{:0.4f}" {:0.3f} deg'.format(sclib[target][band][vislist[0]][selfcal_plan[target][band]['solints'][i]]['Beam_major_pre'],sclib[target][band][vislist[0]][selfcal_plan[target][band]['solints'][i]]['Beam_minor_pre'],sclib[target][band][vislist[0]][selfcal_plan[target][band]['solints'][i]]['Beam_PA_pre'])+'

\n') + htmlOutSolint.writelines('Post SC SNR: {:0.3f}'.format(sclib[target][band][representative_vislist[0]][selfcal_plan[target][band]['solints'][i]]['SNR_post'])+'
Pre SC SNR: {:0.3f}'.format(sclib[target][band][representative_vislist[0]][selfcal_plan[target][band]['solints'][i]]['SNR_pre'])+'

\n') + htmlOutSolint.writelines('Post SC RMS: {:0.7f}'.format(sclib[target][band][representative_vislist[0]][selfcal_plan[target][band]['solints'][i]]['RMS_post'])+' Jy/beam
Pre SC RMS: {:0.7f}'.format(sclib[target][band][representative_vislist[0]][selfcal_plan[target][band]['solints'][i]]['RMS_pre'])+' Jy/beam
\n') + htmlOutSolint.writelines('Post Beam: {:0.4f}"x{:0.4f}" {:0.3f} deg'.format(sclib[target][band][representative_vislist[0]][selfcal_plan[target][band]['solints'][i]]['Beam_major_post'],sclib[target][band][representative_vislist[0]][selfcal_plan[target][band]['solints'][i]]['Beam_minor_post'],sclib[target][band][representative_vislist[0]][selfcal_plan[target][band]['solints'][i]]['Beam_PA_post'])+'
\n') + htmlOutSolint.writelines('Pre Beam: {:0.4f}"x{:0.4f}" {:0.3f} deg'.format(sclib[target][band][representative_vislist[0]][selfcal_plan[target][band]['solints'][i]]['Beam_major_pre'],sclib[target][band][representative_vislist[0]][selfcal_plan[target][band]['solints'][i]]['Beam_minor_pre'],sclib[target][band][representative_vislist[0]][selfcal_plan[target][band]['solints'][i]]['Beam_PA_pre'])+'

\n') if 'inf_EB' in selfcal_plan[target][band]['solints'][i]: htmlOutSolint.writelines('

Phase vs. Frequency Plots:

\n') else: htmlOutSolint.writelines('

Phase vs. Time Plots:

\n') - for vis in vislist: + for vis in representative_vislist: htmlOutSolint.writelines('

MS: '+vis+'

\n') - if 'gaintable' not in sclib[target][band][vis][selfcal_plan[target][band]['solints'][i]]: + if selfcal_plan[target][band]['solints'][i] not in sclib[target][band][vis] or 'gaintable' not in sclib[target][band][vis][selfcal_plan[target][band]['solints'][i]]: htmlOutSolint.writelines('No gaintable available

') continue ant_list=get_ant_list(vis)