diff --git a/auto_selfcal/__main__.py b/auto_selfcal/__main__.py
index 42ce51d5..cc195108 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 a475e76c..3eebbcd3 100644
--- a/auto_selfcal/auto_selfcal.py
+++ b/auto_selfcal/auto_selfcal.py
@@ -8,6 +8,7 @@
import sys
import pickle
import pprint
+import copy
from .selfcal_helpers import *
from .run_selfcal import run_selfcal
@@ -15,6 +16,7 @@
from .weblog_creation import *
from .prepare_selfcal import prepare_selfcal, set_clean_thresholds, plan_selfcal_per_solint
from .original_ms_helpers import applycal_to_orig_MSes, uvcontsub_orig_MSes
+from .prepare_cocal import prepare_cocal
try:
from .alignment_helpers import align_measurement_sets
except:
@@ -36,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,
@@ -449,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)
@@ -630,6 +633,7 @@ def auto_selfcal(
with open('selfcal_library.pickle', 'wb') as handle:
pickle.dump(selfcal_library, handle, protocol=pickle.HIGHEST_PROTOCOL)
+
import json
class NpEncoder(json.JSONEncoder):
@@ -709,7 +713,8 @@ def default(self, obj):
with open('selfcal_library.pickle', 'wb') as handle:
pickle.dump(selfcal_library, handle, protocol=pickle.HIGHEST_PROTOCOL)
-
+ with open('selfcal_plan.pickle', 'wb') as handle:
+ pickle.dump(selfcal_plan, handle, protocol=pickle.HIGHEST_PROTOCOL)
##
## Begin Self-cal loops
@@ -729,109 +734,7 @@ def default(self, obj):
if allow_cocal:
- ##
- ## Save the flags following the main iteration of self-calibration since we will need to revert to the beginning for the fallback mode.
- ##
- # PS: I don't need this anymore?
- for vis in selfcal_library[target][band]['vislist']:
- if not os.path.exists(vis+'.flagversions/flags.fb_selfcal_starting_flags'):
- flagmanager(vis=vis,mode='save',versionname='fb_selfcal_starting_flags')
- else:
- flagmanager(vis=vis,mode='restore',versionname='fb_selfcal_starting_flags')
-
- ##
- ## For sources that self-calibration failed, try to use the inf_EB and the inf solutions from the sources that
- ## were successful.
-
- for target in selfcal_library.keys():
- for band in selfcal_library[target].keys():
- print(target, selfcal_library[target][band]["final_solint"])
-
- inf_EB_fields = {}
- inf_fields = {}
- fallback_fields = {}
- calibrators = {}
- for band in bands:
- # Initialize the lists for this band.
- inf_EB_fields[band] = []
- inf_fields[band] = []
- fallback_fields[band] = []
-
- # Loop through and identify which sources belong where.
- for target in selfcal_library.keys():
- if selfcal_library[target][band]['SC_success'] and 'fb' not in selfcal_library[target][band]['final_solint']:
- inf_EB_fields[band].append(target)
- if selfcal_library[target][band]['final_solint'] != 'inf_EB':
- inf_fields[band].append(target)
- elif 'inf' in selfcal_plan[target][band]['solints']:
- fallback_fields[band].append(target)
- else:
- fallback_fields[band].append(target)
-
- # Update the relevant lists if we are going to do a fallback mode.
- if len(fallback_fields[band]) > 0:
- 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]['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]]
- applycal_mode[band][target] += [applycal_mode[band][target][0], applycal_mode[band][target][1], applycal_mode[band][target][1], applycal_mode[band][target][1]]
- calibrators[band] = [inf_EB_fields[band], inf_fields[band], inf_fields[band], inf_fields[band]]
- selfcal_library[target][band]["nsigma"] = np.concatenate((selfcal_library[target][band]["nsigma"],[selfcal_library[target][band]["nsigma"][0], \
- selfcal_library[target][band]["nsigma"][1], selfcal_library[target][band]["nsigma"][1], selfcal_library[target][band]["nsigma"][1]]))
-
- print(inf_EB_fields)
- print(inf_fields)
- print(fallback_fields)
-
- ##
- ## Reset the inf_EB informational dictionaries.
- ##
-
- for target in selfcal_library:
- for band in solint_snr[target].keys():
- # If the target had a successful inf_EB solution, no need to reset.
- if target in inf_EB_fields[band]:
- continue
-
- for vis in selfcal_library[target][band]['vislist']:
- 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'
- selfcal_plan[target][band][vis]['inf_EB_gaintype']=inf_EB_gaintype #G
- selfcal_plan[target][band][vis]['inf_EB_fallback_mode']='' #'scan'
-
-
- calculate_inf_EB_fb_anyways = True
- preapply_targets_own_inf_EB = False
-
- ## The below sets the calibrations back to what they were prior to starting the fallback mode. It should not be needed
- ## for the final version of the codue, but is used for testing.
-
-
- for target in selfcal_library:
- sani_target=sanitize_string(target)
- for band in selfcal_library[target].keys():
- if target not in fallback_fields[band]:
- continue
- if 'gaintable_final' in selfcal_library[target][band][vislist[0]]:
- print('****************Reapplying previous solint solutions*************')
- for vis in selfcal_library[target][band]['vislist']:
- print('****************Applying '+str(selfcal_library[target][band][vis]['gaintable_final'])+' to '+target+' '+band+'*************')
- ## NOTE: should this be selfcal_starting_flags instead of fb_selfcal_starting_flags ???
- flagmanager(vis=vis,mode='delete',versionname='fb_selfcal_starting_flags_'+sani_target)
- applycal(vis=vis,\
- gaintable=selfcal_library[target][band][vis]['gaintable_final'],\
- interp=selfcal_library[target][band][vis]['applycal_interpolate_final'],\
- calwt=True,spwmap=selfcal_library[target][band][vis]['spwmap_final'],\
- applymode=selfcal_library[target][band][vis]['applycal_mode_final'],\
- field=target,spw=selfcal_library[target][band][vis]['spws'])
- else:
- print('****************Removing all calibrations for '+target+' '+band+'**************')
- for vis in selfcal_library[target][band]['vislist']:
- flagmanager(vis=vis,mode='delete',versionname='fb_selfcal_starting_flags_'+sani_target)
- clearcal(vis=vis,field=target,spw=selfcal_library[target][band][vis]['spws'])
- ## END
-
-
+ fallback_fields, calibrators = prepare_cocal(selfcal_library, selfcal_plan, inf_EB_gaincal_combine, inf_EB_gaintype)
##
## Begin fallback self-cal loops
##
@@ -845,8 +748,7 @@ def default(self, obj):
inf_EB_gaincal_combine=inf_EB_gaincal_combine, inf_EB_gaintype=inf_EB_gaintype, unflag_only_lbants=unflag_only_lbants, \
unflag_only_lbants_onlyap=unflag_only_lbants_onlyap, calonly_max_flagged=calonly_max_flagged, \
second_iter_solmode=second_iter_solmode, unflag_fb_to_prev_solint=unflag_fb_to_prev_solint, rerank_refants=rerank_refants, \
- mode="cocal", calibrators=calibrators, calculate_inf_EB_fb_anyways=calculate_inf_EB_fb_anyways, \
- preapply_targets_own_inf_EB=preapply_targets_own_inf_EB, gaincalibrator_dict=gaincalibrator_dict, allow_gain_interpolation=True)
+ mode="cocal", calibrators=calibrators, gaincalibrator_dict=gaincalibrator_dict, allow_gain_interpolation=True)
if debug:
print(json.dumps(selfcal_library, indent=4, cls=NpEncoder))
diff --git a/auto_selfcal/gaincal_wrapper.py b/auto_selfcal/gaincal_wrapper.py
index dab9e4cf..3cf72b67 100644
--- a/auto_selfcal/gaincal_wrapper.py
+++ b/auto_selfcal/gaincal_wrapper.py
@@ -1,11 +1,12 @@
import numpy as np
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,
- guess_scan_combine=False, do_fallback_calonly=False):
+ guess_scan_combine=False, do_fallback_calonly=False, preapply_targets_own_inf_EB=False):
"""
This function runs gaincal for a given target, band, and solint, and updates the selfcal_library and selfcal_plan dictionaries with the results.
It also handles the pre-application of inf_EB solutions if necessary.
@@ -26,11 +27,14 @@ 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.
- for mode in selfcal_plan[vis]['solint_settings'][solint]['modes_to_attempt']:
- selfcal_plan[vis]['solint_settings'][solint]['gaincal_return_dict'][mode] = []
+ for gc_mode in selfcal_plan[vis]['solint_settings'][solint]['modes_to_attempt']:
+ selfcal_plan[vis]['solint_settings'][solint]['gaincal_return_dict'][gc_mode] = []
+
##
## Set gaincal parameters depending on which iteration and whether to use combine=spw for inf_EB or not
## Defaults should assume combine='scan' and gaintpe='G' will fallback to combine='scan,spw' if too much flagging
@@ -39,58 +43,40 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so
current_solint_index=selfcal_plan['solints'].index(solint)
selfcal_plan[vis]['solint_settings'][solint]['computed_gaintable'] = {}
if selfcal_plan['solmode'][iteration] == 'p':
- if mode == "cocal":
- if 'inf_EB' in selfcal_library[vis]:
- #gaincal_preapply_gaintable[vis]=[sani_target+'_'+vis+'_'+band+'_inf_EB_0_p.g']
- if calculate_inf_EB_fb_anyways or not selfcal_library[vis]["inf_EB"]["Pass"]:
- previous_solint = "inf_EB_fb"
- else:
- previous_solint = "inf_EB"
- else:
- #gaincal_preapply_gaintable[vis]=[sani_target+'_'+vis+'_'+band+'_inf_EB_fb_'+str(iteration-1)+'_p.g']
- previous_solint = "inf_EB_fb"
- else:
- # not entirely sure why this if/else is necessary might not be with new selfcal plan
- if selfcal_plan['solmode'][iteration]=='p':
- previous_solint = "inf_EB"
- else:
- previous_solint = selfcal_library['final_phase_solint']
gaincal_spwmap=[]
gaincal_preapply_gaintable=[]
gaincal_interpolate=[]
applycal_spwmap=[]
applycal_interpolate=[]
applycal_gaintable=[]
- for j in range(current_solint_index):
- if selfcal_plan[vis]['solint_settings'][selfcal_plan['solints'][j]]['preapply_this_gaintable']:
- gaincal_preapply_gaintable.append(selfcal_plan[vis]['solint_settings'][selfcal_plan['solints'][j]]['accepted_gaintable'])
- gaincal_spwmap.append(selfcal_plan[vis]['solint_settings'][selfcal_plan['solints'][j]]['applycal_spwmap'])
- gaincal_interpolate.append(selfcal_plan[vis]['solint_settings'][selfcal_plan['solints'][j]]['applycal_interpolate'])
- applycal_spwmap.append(selfcal_plan[vis]['solint_settings'][selfcal_plan['solints'][j]]['applycal_spwmap'])
- applycal_interpolate.append(selfcal_plan[vis]['solint_settings'][selfcal_plan['solints'][j]]['applycal_interpolate'])
- applycal_gaintable.append(selfcal_plan[vis]['solint_settings'][selfcal_plan['solints'][j]]['accepted_gaintable'])
+
+ if mode == "cocal":
+ include_solints = selfcal_plan[vis]['solint_settings'][solint]['preapply_solints']
+ else:
+ include_solints = [selfcal_plan['solints'][j] for j in range(current_solint_index) if
+ selfcal_plan[vis]['solint_settings'][selfcal_plan['solints'][j]]['preapply_this_gaintable']]
+
+ for incl_solint in include_solints:
+ gaincal_preapply_gaintable.append(selfcal_plan[vis]['solint_settings'][incl_solint]['accepted_gaintable'])
+ gaincal_spwmap.append(selfcal_plan[vis]['solint_settings'][incl_solint]['applycal_spwmap'])
+ gaincal_interpolate.append(selfcal_plan[vis]['solint_settings'][incl_solint]['applycal_interpolate'])
+
+ # If we manually specify that applycal should use a different solint reference than gaincal (cocal), change
+ # to that solint here.
+ if "applycal_solint" in selfcal_plan[vis]['solint_settings'][incl_solint]:
+ incl_solint = selfcal_plan[vis]['solint_settings'][incl_solint]["applycal_solint"]
+
+ applycal_spwmap.append(selfcal_plan[vis]['solint_settings'][incl_solint]['applycal_spwmap'])
+ applycal_interpolate.append(selfcal_plan[vis]['solint_settings'][incl_solint]['applycal_interpolate'])
+ if solint == "inf_fb3":
+ applycal_gaintable.append(selfcal_plan[vis]['solint_settings']["inf_EB_fb"]['accepted_gaintable'])
+ else:
+ applycal_gaintable.append(selfcal_plan[vis]['solint_settings'][incl_solint]['accepted_gaintable'])
if solint != 'inf_EB':
gaincal_solmode = "" if not do_fallback_calonly or second_iter_solmode == "GSPLINE" else second_iter_solmode
- """
- if 'spw' in selfcal_plan[vis]['inf_EB_gaincal_combine']:
- applycal_spwmap[vis]=[selfcal_library[vis]['spwmap'],selfcal_library[vis]['spwmap']]
- gaincal_spwmap[vis]=[selfcal_library[vis]['spwmap']]
- elif selfcal_plan[vis]['inf_EB_fallback_mode']=='spwmap':
- applycal_spwmap[vis]=selfcal_library[vis]['inf_EB']['spwmap'] + [selfcal_library[vis]['spwmap']]
- gaincal_spwmap[vis]=selfcal_library[vis]['inf_EB']['spwmap']
- else:
- applycal_spwmap[vis]=[[],selfcal_library[vis]['spwmap']]
- gaincal_spwmap[vis]=[]
- """
- # Revert back to applying the inf_EB solution if calculate_inf_EB_fb_anyways, i.e. we just use the inf_EB_fb solution
- # for gaincal.
- if mode == "cocal":
- if selfcal_library['final_solint'] == 'inf_EB' and calculate_inf_EB_fb_anyways:
- previous_solint = "inf_EB"
-
fallback=''
if selfcal_plan['solmode'][iteration] == 'ap':
solnorm=True
@@ -101,266 +87,118 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so
# Check which targets are acceptable to use as calibrators.
targets = calibrators[band][iteration - len(selfcal_plan['solints'])]
- include_targets, include_scans = triage_calibrators(vis, target, targets)
+ include_targets, include_scans = triage_calibrators(vis, target, band, targets)
else:
#include_targets = str(selfcal_library['sub-fields-fid_map'][vis][0])
include_targets = selfcal_library['bands_for_targets'][vis]['field_str']
include_scans = ""
+ print('Include targets: ', include_targets)
- if solint == "scan_inf":
- if len(gaincalibrator_dict[vis]) > 0:
- print("Determining scan_inf from calibrator scans in full MS")
- scans = []
- intents = []
- times = []
- for t in gaincalibrator_dict[vis].keys():
- scans += [gaincalibrator_dict[vis][t]["scans"]]
- intents += [np.repeat(gaincalibrator_dict[vis][t]["intent"],gaincalibrator_dict[vis][t]["scans"].size)]
- times += [gaincalibrator_dict[vis][t]["times"]]
-
- times = np.concatenate(times)
- order = np.argsort(times)
- times = times[order]
-
- scans = np.concatenate(scans)[order]
- intents = np.concatenate(intents)[order]
-
- is_gaincalibrator = intents == "phase"
- scans = scans[is_gaincalibrator]
-
- msmd.open(vis)
- if selfcal_library['telescope'] == 'ALMA' or selfcal_library['telescope'] == 'ACA':
- scan_ids_for_target = msmd.scansforfield(target)
- elif 'VLA' in selfcal_library['telescope']:
- scan_ids_for_target=np.array([],dtype=int)
- for fid in selfcal_library['sub-fields']:
- if fid in selfcal_library['sub-fields-fid_map'][vis].keys():
- field_id=selfcal_library['sub-fields-fid_map'][vis][fid]
- scan_ids_for_target=np.append(scan_ids_for_target,msmd.scansforfield(field_id))
-
- scan_ids_for_target.sort() # sort scans since they will be out of order
- include_scans = []
- for iscan in range(scans.size-1):
- scan_group = np.intersect1d(scan_ids_for_target,
- np.array(list(range(scans[iscan]+1, scans[iscan+1])))).astype(str)
- if scan_group.size > 0:
- include_scans.append(",".join(scan_group))
- # PIPE-2741: if there are any scans before the first gain calibrator scan or after the last gain calibrator scan, catch them.
- if scans.size > 0:
- extra_scans = scan_ids_for_target[scan_ids_for_target > max(scans)]
- if extra_scans.size > 0:
- include_scans.append(','.join(extra_scans.astype(str)))
- extra_scans = scan_ids_for_target[scan_ids_for_target < min(scans)]
- if extra_scans.size > 0:
- include_scans.append(','.join(extra_scans.astype(str)))
- msmd.close()
- elif guess_scan_combine:
- print("Determining scan_inf from guessing where the calibrator scans were")
- msmd.open(vis)
- include_scans = []
-
- #to guess at scan_inf combination for VLA look for breaks in the consecutive
- #scan numbers and assume that the break is due to a calibrator scan
- #Fetch scans for scan inf by collecting the field ids and running msmd.scansforfield
- if 'VLA' in selfcal_library['telescope']:
- scans=np.array([],dtype=int)
- for fid in selfcal_library['sub-fields']:
- if fid in selfcal_library['sub-fields-fid_map'][vis].keys():
- field_id=selfcal_library['sub-fields-fid_map'][vis][fid]
- scans=np.append(scans,msmd.scansforfield(field_id))
-
- scans.sort() # sort scans since they will be out of order
- scan_group=''
- for iscan in range(scans.size):
- if len(include_scans) > 0:
- if str(scans[iscan]) in include_scans[-1]:
- continue
- if scan_group == '':
- scan_group = str(int(scans[iscan]))
-
- if iscan < scans.size-1:
- if scans[iscan+1] == scans[iscan]+1:
- scan_group += ","+str(int(scans[iscan+1]))
- else:
- include_scans.append(scan_group)
- scan_group=''
- else: #write out the last scan group to include_scans
- if scan_group != '':
- include_scans.append(scan_group)
- scan_group=''
-
- #guess scan_inf combination by getting all the scans for targets and do a simple grouping
- if selfcal_library['telescope'] == 'ALMA' or selfcal_library['telescope'] == 'ACA':
- scans = msmd.scansforfield(target)
-
- for iscan in range(scans.size):
- if len(include_scans) > 0:
- if str(scans[iscan]) in include_scans[-1]:
- continue
-
- scan_group = str(scans[iscan])
-
- if iscan < scans.size-1:
- if msmd.fieldsforscan(scans[iscan+1]).size < msmd.fieldsforscan(scans[iscan]).size/3:
- scan_group += ","+str(scans[iscan+1])
-
- include_scans.append(scan_group)
-
- msmd.close()
- else:
- print("Not guessing where calibration scans are and justincluding all scans")
- msmd.open(vis)
- if selfcal_library['telescope'] == 'ALMA' or selfcal_library['telescope'] == 'ACA':
- include_scans = [str(scan) for scan in msmd.scansforfield(target)]
- elif 'VLA' in selfcal_library['telescope']:
- scans=np.array([],dtype=int)
- for fid in selfcal_library['sub-fields']:
- if fid in selfcal_library['sub-fields-fid_map'][vis].keys():
- field_id=selfcal_library['sub-fields-fid_map'][vis][fid]
- scans=np.append(scans,msmd.scansforfield(field_id))
-
- scans.sort() # sort scans since they will be out of order
- include_scans = [str(scan) for scan in scans]
- msmd.close()
+ 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]
- if mode == "cocal" and preapply_targets_own_inf_EB and "inf_fb" in solint and "inf" in selfcal_plan['solints']:
- ##
- ## If we want to pre-apply inf_EB solution from each calibrator to itself, all we do is combine all of thier
- ## individual inf tables, as these were pre-calculated in that way.
- ##
- destination_table = sani_target+'_'+vis+'_'+band+'_'+solint+'_'+str(iteration)+'_'+selfcal_plan['solmode'][iteration]+'.g'
- for t in include_targets.split(","):
- if os.path.exists(sanitize_string(t)+'_'+vis+'_'+band+'_'+solint.replace('_fb1','').replace('_fb2','').replace('_fb3','')+'_'+str(1)+'_'+selfcal_plan['solmode'][iteration]+\
- '.pre-pass.g'):
- table_name = sanitize_string(t)+'_'+vis+'_'+band+'_'+solint.replace('_fb1','').replace('_fb2','').replace('_fb3','')+'_'+str(1)+'_'+selfcal_plan['solmode'][iteration]+\
- '.pre-pass.g'
+
+ # Fields that don't have any mask in the primary beam should be removed from consideration, as their models are likely bad.
+ gaincal_solmode=""
+ if selfcal_library['obstype'] == 'mosaic':
+ msmd.open(vis)
+ include_targets = []
+ remove = []
+ 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[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[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)
+
+ if len(scan_targets) > 0:
+ include_targets.append(','.join(np.array(scan_targets).astype(str)))
else:
- table_name = sanitize_string(t)+'_'+vis+'_'+band+'_'+solint.replace('_fb1','').replace('_fb2','').replace('_fb3','')+'_'+str(1)+'_'+selfcal_plan['solmode'][iteration]+'.g'
- #t_final_solint = selfcal_library[t][band]["final_phase_solint"]
- #t_iteration = selfcal_library[t][band][vislist[0]][t_final_solint]["iteration"]
- #table_name = sanitize_string(t)+'_'+vis+'_'+band+'_'+t_final_solint+'_'+str(t_iteration)+'_'+solmode[band][iteration]+'.g'
-
- rerefant(vis, table_name, caltable="tmp0.g", refantmode="strict", refant=selfcal_library[vis]['refant'])
+ remove.append(incl_scan)
- tb.open("tmp0.g")
- if not os.path.exists("tmp1.g"):
- tb.copy("tmp1.g", deep=True)
- else:
- tb.copyrows("tmp1.g")
- tb.close()
+ for incl_scan in remove:
+ include_scans.remove(incl_scan)
- os.system("rm -rf tmp0.g")
+ msmd.close()
+ elif solint == "inf_fb3":
+ include_targets = include_targets.split(',')
+ include_scans = include_scans * len(include_targets)
+ else:
+ include_targets = [include_targets]
+
+ selfcal_library[vis][solint]["include_scans"] = include_scans
+ selfcal_library[vis][solint]["include_targets"] = include_targets
+ print(solint,'Include scans: ', include_scans)
+ print(solint,'Include targets: ', include_targets)
+ print(solint,'Modes to attempt: ',selfcal_plan[vis]['solint_settings'][solint]['modes_to_attempt'])
+ for incl_scans, incl_targets in zip(include_scans, include_targets):
+ if solint == "inf_fb3":
+ gaincal_preapply_gaintable = [selfcal_plan[vis]['solint_settings']["inf_fb3"]["preapply_gaintable_dict"][incl_targets]]
+
+ for gc_mode in selfcal_plan[vis]['solint_settings'][solint]['modes_to_attempt']:
+ print(incl_targets)
+ print(incl_scans)
+ print(vis,solint,gc_mode)
+ print(selfcal_plan[vis]['solint_settings'][solint]['gaincal_combine'])
+ gaincal_combine=selfcal_plan[vis]['solint_settings'][solint]['gaincal_combine'][gc_mode]
+ filename_append=selfcal_plan[vis]['solint_settings'][solint]['filename_append'][gc_mode]
- tb.open("tmp1.g")
- subt = tb.query("OBSERVATION_ID==0", sortlist="TIME,ANTENNA1")
- copyt = subt.copy(destination_table, deep=True)
- tb.close()
- subt.close()
- copyt.close()
+ if 'spw' in gaincal_combine:
+ if selfcal_library['spws_set'][vis].ndim == 1:
+ nspw_sets=1
+ else:
+ nspw_sets=selfcal_library['spws_set'][vis].shape[0]
+ else: #only necessary to loop over gain cal when in inf_EB to avoid inf_EB solving for all spws
+ nspw_sets=1
- os.system("rm -rf tmp1.g")
+ for i in range(nspw_sets): # run gaincal on each spw set to handle spectral scans
+ if 'spw' in gaincal_combine:
+ if nspw_sets == 1 and selfcal_library['spws_set'][vis].ndim == 1:
+ spwselect=','.join(str(spw) for spw in selfcal_library['spws_set'][vis].tolist())
+ else:
+ spwselect=','.join(str(spw) for spw in selfcal_library['spws_set'][vis][i].tolist())
+ else:
+ spwselect=selfcal_library[vis]['spws']
- # Remove all of the scans that failed the triage above.
- tb.open(destination_table, nomodify=False)
- scans = tb.getcol("SCAN_NUMBER")
+ gaintable_name=sani_target+'_'+vis+'_'+band+'_'+solint+'_'+str(iteration)+'_'+selfcal_plan['solmode'][iteration]+'_'+filename_append+'.g'
+ print('prior to gaincal',gaintable_name, gc_mode)
- bad_scans = np.repeat(True, scans.size)
- for scan in include_scans[0].split(","):
- bad_scans[scans == int(scan)] = False
+ 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=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',\
+ append=os.path.exists(sani_target+'_'+vis+'_'+band+'_'+solint+'_'+str(iteration)+'_'+selfcal_plan['solmode'][iteration]+'_'+filename_append+'.g'))
+ selfcal_plan[vis]['solint_settings'][solint]['gaincal_return_dict'][gc_mode].append(gcdict.copy())
+ else:
+ for baseband in selfcal_library[vis]['baseband'].keys():
+ 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=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',\
+ append=os.path.exists(sani_target+'_'+vis+'_'+band+'_'+solint+'_'+str(iteration)+'_'+selfcal_plan['solmode'][iteration]+'_'+filename_append+'.g'))
+ selfcal_plan[vis]['solint_settings'][solint]['gaincal_return_dict'][gc_mode].append(gcdict.copy())
- tb.removerows(rownrs=np.where(bad_scans)[0])
- tb.flush()
- tb.close()
- else:
- # Fields that don't have any mask in the primary beam should be removed from consideration, as their models are likely bad.
- gaincal_solmode=""
- if selfcal_library['obstype'] == 'mosaic':
- msmd.open(vis)
- include_targets = []
- remove = []
- 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(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()))]):
- # Note: because of the msmd above getting actual fids from the MS, we just need to append fid below.
- scan_targets.append(fid)
-
- if len(scan_targets) > 0:
- include_targets.append(','.join(np.array(scan_targets).astype(str)))
- else:
- remove.append(incl_scan)
+ selfcal_plan[vis]['solint_settings'][solint]['computed_gaintable'][gc_mode] = gaintable_name
- for incl_scan in remove:
- include_scans.remove(incl_scan)
+ # restricted gaincal table comparisons to only inf_EB prior to changes
+ # commenting because we want to do comparisons for other solints as well
+ #if 'inf_EB' not in solint:
+ # break
- msmd.close()
- else:
- include_targets = [include_targets]
-
- selfcal_library[vis][solint]["include_scans"] = include_scans
- selfcal_library[vis][solint]["include_targets"] = include_targets
- print(solint,'Include scans: ', include_scans)
- print(solint,'Include targets: ', include_targets)
- print(solint,'Modes to attempt: ',selfcal_plan[vis]['solint_settings'][solint]['modes_to_attempt'])
- for incl_scans, incl_targets in zip(include_scans, include_targets):
- for mode in selfcal_plan[vis]['solint_settings'][solint]['modes_to_attempt']:
- print(vis,solint,mode)
- print(selfcal_plan[vis]['solint_settings'][solint]['gaincal_combine'])
- gaincal_combine=selfcal_plan[vis]['solint_settings'][solint]['gaincal_combine'][mode]
- filename_append=selfcal_plan[vis]['solint_settings'][solint]['filename_append'][mode]
-
- if 'spw' in gaincal_combine:
- if selfcal_library['spws_set'][vis].ndim == 1:
- nspw_sets=1
- else:
- nspw_sets=selfcal_library['spws_set'][vis].shape[0]
- else: #only necessary to loop over gain cal when in inf_EB to avoid inf_EB solving for all spws
- nspw_sets=1
- for i in range(nspw_sets): # run gaincal on each spw set to handle spectral scans
- if 'spw' in gaincal_combine:
- if nspw_sets == 1 and selfcal_library['spws_set'][vis].ndim == 1:
- spwselect=','.join(str(spw) for spw in selfcal_library['spws_set'][vis].tolist())
- else:
- 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, mode)
- if mode != 'per_bb':
- gcdict=call_gaincal(vis=vis, caltable=gaintable_name, gaintype=selfcal_plan[vis]['solint_settings'][solint]['gaincal_gaintype'][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',''),\
- 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',\
- append=os.path.exists(sani_target+'_'+vis+'_'+band+'_'+solint+'_'+str(iteration)+'_'+selfcal_plan['solmode'][iteration]+'_'+filename_append+'.g'))
- selfcal_plan[vis]['solint_settings'][solint]['gaincal_return_dict'][mode].append(gcdict.copy())
- else:
- for baseband in selfcal_library[vis]['baseband'].keys():
- 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'][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',''),\
- 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',\
- append=os.path.exists(sani_target+'_'+vis+'_'+band+'_'+solint+'_'+str(iteration)+'_'+selfcal_plan['solmode'][iteration]+'_'+filename_append+'.g'))
- selfcal_plan[vis]['solint_settings'][solint]['gaincal_return_dict'][mode].append(gcdict.copy())
-
- selfcal_plan[vis]['solint_settings'][solint]['computed_gaintable'][mode] = gaintable_name
-
- # restricted gaincal table comparisons to only inf_EB prior to changes
- # commenting because we want to do comparisons for other solints as well
- #if 'inf_EB' not in solint:
- # break
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 = \
@@ -409,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]+''
@@ -419,6 +257,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so
selfcal_library[fid][vis][solint]['applycal_mode']=selfcal_plan['applycal_mode'][iteration]+''
selfcal_library[fid][vis][solint]['applycal_interpolate']=applycal_interpolate.copy()
selfcal_library[fid][vis][solint]['solmode']=selfcal_plan['solmode'][iteration]+''
+
selfcal_plan[vis]['solint_settings'][solint]['accepted_gaintable']=sani_target+'_'+vis+'_'+band+'_'+solint+'_'+str(iteration)+'_'+selfcal_plan['solmode'][iteration]+'_'+preferred_mode+'.g'
if 'combinespw' not in preferred_mode: # preapply all non spwcombine gain tables
selfcal_plan[vis]['solint_settings'][solint]['preapply_this_gaintable']=True
@@ -427,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=[]
@@ -435,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'])
@@ -451,9 +292,9 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so
else:
solnorm=False
- for mode in selfcal_plan[vis]['solint_settings'][solint]['modes_to_attempt']:
- gaincal_combine=selfcal_plan[vis]['solint_settings'][solint]['gaincal_combine'][mode]
- filename_append=selfcal_plan[vis]['solint_settings'][solint]['filename_append'][mode]
+ for gc_mode in selfcal_plan[vis]['solint_settings'][solint]['modes_to_attempt']:
+ gaincal_combine=selfcal_plan[vis]['solint_settings'][solint]['gaincal_combine'][gc_mode]
+ filename_append=selfcal_plan[vis]['solint_settings'][solint]['filename_append'][gc_mode]
if 'spw' in gaincal_combine:
if selfcal_library['spws_set'][vis].ndim == 1:
@@ -470,29 +311,30 @@ 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 mode != 'per_bb':
- gcdict=call_gaincal(vis=vis, caltable=gaintable_name, gaintype=selfcal_plan[vis]['solint_settings'][solint]['gaincal_gaintype'][mode], spw=spwselect,
+ 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',\
append=os.path.exists(gaintable_name))
- selfcal_plan[vis]['solint_settings'][solint]['gaincal_return_dict'][mode].append(gcdict.copy())
+ selfcal_plan[vis]['solint_settings'][solint]['gaincal_return_dict'][gc_mode].append(gcdict.copy())
else:
for baseband in selfcal_library[vis]['baseband'].keys():
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'][mode], spw=spwselect_bb,
+ 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',\
append=os.path.exists(gaintable_name))
- selfcal_plan[vis]['solint_settings'][solint]['gaincal_return_dict'][mode].append(gcdict.copy())
+ selfcal_plan[vis]['solint_settings'][solint]['gaincal_return_dict'][gc_mode].append(gcdict.copy())
- selfcal_plan[vis]['solint_settings'][solint]['computed_gaintable'][mode] = gaintable_name
+ selfcal_plan[vis]['solint_settings'][solint]['computed_gaintable'][gc_mode] = gaintable_name
gaintable_prefix='temp_'
if len(selfcal_plan[vis]['solint_settings'][solint]['modes_to_attempt']) > 1:
@@ -501,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=[]
@@ -531,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]+''
@@ -548,18 +390,18 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, so
filename_append=selfcal_plan[vis]['solint_settings'][solint]['filename_append'][preferred_mode]
- for mode in selfcal_plan[vis]['solint_settings'][solint]['modes_to_attempt']:
+ for gc_mode in selfcal_plan[vis]['solint_settings'][solint]['modes_to_attempt']:
##### send chosen subtable to this routine for final copying to the gain table we want.
- tb.open('temp_'+solint+'_'+str(iteration)+'_'+selfcal_plan['solmode'][iteration]+'_'+mode+'.g')
+ tb.open('temp_'+solint+'_'+str(iteration)+'_'+selfcal_plan['solmode'][iteration]+'_'+gc_mode+'.g')
subt = tb.query("OBSERVATION_ID==0", sortlist="TIME,ANTENNA1")
tb.close()
- subt.copy(sani_target+'_'+vis+'_'+band+'_'+solint+'_'+str(iteration)+'_'+selfcal_plan['solmode'][iteration]+'_'+mode+'.g', deep=True)
+ subt.copy(sani_target+'_'+vis+'_'+band+'_'+solint+'_'+str(iteration)+'_'+selfcal_plan['solmode'][iteration]+'_'+gc_mode+'.g', deep=True)
subt.close()
# commented so that we keep all the gain tables around
# once well-tested, we might remove the tables for the non-chosen modes to avoid generating too many useless files.
- os.system('rm -rf '+'temp_'+solint+'_'+str(iteration)+'_'+selfcal_plan['solmode'][iteration]+'_'+mode+'.g')
+ os.system('rm -rf '+'temp_'+solint+'_'+str(iteration)+'_'+selfcal_plan['solmode'][iteration]+'_'+gc_mode+'.g')
if rerank_refants:
selfcal_library[vis]["refant"] = rank_refants(vis, selfcal_library['telescope'], caltable=sani_target+'_'+vis+'_'+band+'_'+solint+'_'+str(iteration)+'_'+selfcal_plan['solmode'][iteration]+'_'+selfcal_library[vis][solint]['final_mode']+'.g')
@@ -606,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.
@@ -650,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()
@@ -673,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'
@@ -708,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 20b0a771..a4690934 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 93543f2a..06715f21 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,147 @@ 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
+
+def scan_inf_scan_combine(selfcal_library, vis, target, gaincalibrator_dict, guess_scan_combine=False):
+ if len(gaincalibrator_dict[vis]) > 0:
+ print("Determining scan_inf from calibrator scans in full MS")
+ scans = []
+ intents = []
+ times = []
+ for t in gaincalibrator_dict[vis].keys():
+ scans += [gaincalibrator_dict[vis][t]["scans"]]
+ intents += [np.repeat(gaincalibrator_dict[vis][t]["intent"],gaincalibrator_dict[vis][t]["scans"].size)]
+ times += [gaincalibrator_dict[vis][t]["times"]]
+
+ times = np.concatenate(times)
+ order = np.argsort(times)
+ times = times[order]
+
+ scans = np.concatenate(scans)[order]
+ intents = np.concatenate(intents)[order]
+
+ is_gaincalibrator = intents == "phase"
+ scans = scans[is_gaincalibrator]
+
+ msmd.open(vis)
+ if selfcal_library['telescope'] == 'ALMA' or selfcal_library['telescope'] == 'ACA':
+ scan_ids_for_target = msmd.scansforfield(target)
+ elif 'VLA' in selfcal_library['telescope']:
+ scan_ids_for_target=np.array([],dtype=int)
+ for fid in selfcal_library['sub-fields']:
+ if fid in selfcal_library['sub-fields-fid_map'][vis].keys():
+ field_id=selfcal_library['sub-fields-fid_map'][vis][fid]
+ scan_ids_for_target=np.append(scan_ids_for_target,msmd.scansforfield(field_id))
+
+ scan_ids_for_target.sort() # sort scans since they will be out of order
+ include_scans = []
+ for iscan in range(scans.size-1):
+ scan_group = np.intersect1d(scan_ids_for_target,
+ np.array(list(range(scans[iscan]+1, scans[iscan+1])))).astype(str)
+ if scan_group.size > 0:
+ include_scans.append(",".join(scan_group))
+ # PIPE-2741: if there are any scans before the first gain calibrator scan or after the last gain calibrator scan, catch them.
+ if scans.size > 0:
+ extra_scans = scan_ids_for_target[scan_ids_for_target > max(scans)]
+ if extra_scans.size > 0:
+ include_scans.append(','.join(extra_scans.astype(str)))
+ extra_scans = scan_ids_for_target[scan_ids_for_target < min(scans)]
+ if extra_scans.size > 0:
+ include_scans.append(','.join(extra_scans.astype(str)))
+ msmd.close()
+ elif guess_scan_combine:
+ print("Determining scan_inf from guessing where the calibrator scans were")
+ msmd.open(vis)
+ include_scans = []
+
+ #to guess at scan_inf combination for VLA look for breaks in the consecutive
+ #scan numbers and assume that the break is due to a calibrator scan
+ #Fetch scans for scan inf by collecting the field ids and running msmd.scansforfield
+ if 'VLA' in selfcal_library['telescope']:
+ scans=np.array([],dtype=int)
+ for fid in selfcal_library['sub-fields']:
+ if fid in selfcal_library['sub-fields-fid_map'][vis].keys():
+ field_id=selfcal_library['sub-fields-fid_map'][vis][fid]
+ scans=np.append(scans,msmd.scansforfield(field_id))
+
+ scans.sort() # sort scans since they will be out of order
+ scan_group=''
+ for iscan in range(scans.size):
+ if len(include_scans) > 0:
+ if str(scans[iscan]) in include_scans[-1]:
+ continue
+ if scan_group == '':
+ scan_group = str(int(scans[iscan]))
+
+ if iscan < scans.size-1:
+ if scans[iscan+1] == scans[iscan]+1:
+ scan_group += ","+str(int(scans[iscan+1]))
+ else:
+ include_scans.append(scan_group)
+ scan_group=''
+ else: #write out the last scan group to include_scans
+ if scan_group != '':
+ include_scans.append(scan_group)
+ scan_group=''
+
+ #guess scan_inf combination by getting all the scans for targets and do a simple grouping
+ if selfcal_library['telescope'] == 'ALMA' or selfcal_library['telescope'] == 'ACA':
+ scans = msmd.scansforfield(target)
+
+ for iscan in range(scans.size):
+ if len(include_scans) > 0:
+ if str(scans[iscan]) in include_scans[-1]:
+ continue
+
+ scan_group = str(scans[iscan])
+
+ if iscan < scans.size-1:
+ if msmd.fieldsforscan(scans[iscan+1]).size < msmd.fieldsforscan(scans[iscan]).size/3:
+ scan_group += ","+str(scans[iscan+1])
+
+ include_scans.append(scan_group)
+
+ msmd.close()
+ else:
+ print("Not guessing where calibration scans are and justincluding all scans")
+ msmd.open(vis)
+ if selfcal_library['telescope'] == 'ALMA' or selfcal_library['telescope'] == 'ACA':
+ include_scans = [str(scan) for scan in msmd.scansforfield(target)]
+ elif 'VLA' in selfcal_library['telescope']:
+ scans=np.array([],dtype=int)
+ for fid in selfcal_library['sub-fields']:
+ if fid in selfcal_library['sub-fields-fid_map'][vis].keys():
+ field_id=selfcal_library['sub-fields-fid_map'][vis][fid]
+ scans=np.append(scans,msmd.scansforfield(field_id))
+
+ scans.sort() # sort scans since they will be out of order
+ include_scans = [str(scan) for scan in scans]
+ msmd.close()
+
+ return include_scans
\ No newline at end of file
diff --git a/auto_selfcal/original_ms_helpers.py b/auto_selfcal/original_ms_helpers.py
index f41bb32f..7a56378e 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
new file mode 100644
index 00000000..cf845663
--- /dev/null
+++ b/auto_selfcal/prepare_cocal.py
@@ -0,0 +1,148 @@
+from casatasks import flagmanager, applycal, clearcal
+from .prepare_selfcal import plan_selfcal_per_solint
+import numpy as np
+import copy
+import os
+
+from .selfcal_helpers import sanitize_string
+
+def prepare_cocal(selfcal_library, selfcal_plan, inf_EB_gaincal_combine, inf_EB_gaintype):
+ for target in selfcal_library:
+ for band in selfcal_library[target]:
+ for vis in selfcal_library[target][band]['vislist']:
+ if not os.path.exists(vis+'.flagversions/flags.fb_selfcal_starting_flags'):
+ flagmanager(vis=vis,mode='save',versionname='fb_selfcal_starting_flags')
+ else:
+ flagmanager(vis=vis,mode='restore',versionname='fb_selfcal_starting_flags')
+
+ ##
+ ## For sources that self-calibration failed, try to use the inf_EB and the inf solutions from the sources that
+ ## were successful.
+
+ for target in selfcal_library.keys():
+ for band in selfcal_library[target].keys():
+ print(target, selfcal_library[target][band]["final_solint"])
+
+ inf_EB_fields = {}
+ inf_fields = {}
+ fallback_fields = {}
+ calibrators = {}
+
+ # Collect the potential targets and calibrators
+ for target in selfcal_library:
+ for band in selfcal_library[target]:
+ if band not in inf_EB_fields:
+ # Initialize the lists for this band.
+ inf_EB_fields[band] = []
+ inf_fields[band] = []
+ fallback_fields[band] = []
+
+ if selfcal_library[target][band]['SC_success'] and 'fb' not in selfcal_library[target][band]['final_solint']:
+ inf_EB_fields[band].append(target)
+ if selfcal_library[target][band]['final_solint'] != 'inf_EB' or \
+ (selfcal_library[target][band]['final_solint'] == 'inf_EB' and
+ 'inf' not in selfcal_plan[target][band]['solints']):
+ inf_fields[band].append(target)
+ else:
+ fallback_fields[band].append(target)
+ else:
+ fallback_fields[band].append(target)
+
+
+ if len(fallback_fields[band]) > 0:
+ for target in selfcal_library:
+ for band in selfcal_library[target]:
+ # 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']:
+ selfcal_plan[target][band][vis]['solint_settings'] = {}
+ for solint_name in ["inf_EB_fb","inf_fb1","inf_fb2","inf_fb3"]:
+ 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,
+ optimize_spw_combine=False,
+ solints=["inf_EB_fb","inf_fb1","inf_fb2","inf_fb3"])
+
+ for target in selfcal_library:
+ for band in selfcal_library[target]:
+ for vis in selfcal_library[target][band]['vislist']:
+ selfcal_plan[target][band][vis]['solint_settings']["inf_EB_fb"]["preapply_solints"] = []
+ selfcal_plan[target][band][vis]['solint_settings']["inf_fb1"]["preapply_solints"] = ["inf_EB_fb"]
+ selfcal_plan[target][band][vis]['solint_settings']["inf_fb2"]["preapply_solints"] = ["inf_EB"]
+ selfcal_plan[target][band][vis]['solint_settings']["inf_fb3"]["preapply_solints"] = ["inf_EB"]
+
+ if selfcal_library[target][band]['SC_success']:
+ selfcal_plan[target][band][vis]['solint_settings']["inf_fb1"]["applycal_solint"] = ["inf_EB"]
+ selfcal_plan[target][band][vis]['solint_settings']["inf_fb2"]["applycal_solint"] = ["inf_EB"]
+ selfcal_plan[target][band][vis]['solint_settings']["inf_fb3"]["applycal_solint"] = ["inf_EB"]
+ else:
+ selfcal_plan[target][band][vis]['solint_settings']["inf_fb1"]["applycal_solint"] = ["inf_EB_fb"]
+ selfcal_plan[target][band][vis]['solint_settings']["inf_fb2"]["applycal_solint"] = ["inf_EB_fb"]
+ selfcal_plan[target][band][vis]['solint_settings']["inf_fb3"]["applycal_solint"] = ["inf_EB_fb"]
+
+ selfcal_plan[target][band][vis]['solint_settings']["inf_fb3"]["preapply_gaintable_dict"] = {}
+ 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]]
+
+ selfcal_library[target][band]["nsigma"] = np.concatenate((selfcal_library[target][band]["nsigma"],[selfcal_library[target][band]["nsigma"][0], \
+ selfcal_library[target][band]["nsigma"][1], selfcal_library[target][band]["nsigma"][1], selfcal_library[target][band]["nsigma"][1]]))
+
+ print(inf_EB_fields)
+ print(inf_fields)
+ print(fallback_fields)
+
+ ##
+ ## Reset the inf_EB informational dictionaries.
+ ##
+
+ for target in selfcal_library:
+ for band in selfcal_library[target].keys():
+ # If the target had a successful inf_EB solution, no need to reset.
+ if target in inf_EB_fields[band]:
+ continue
+
+ for vis in selfcal_library[target][band]['vislist']:
+ 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'
+ selfcal_plan[target][band][vis]['inf_EB_gaintype']=inf_EB_gaintype #G
+ selfcal_plan[target][band][vis]['inf_EB_fallback_mode']='' #'scan'
+
+ ## The below sets the calibrations back to what they were prior to starting the fallback mode. It should not be needed
+ ## for the final version of the codue, but is used for testing.
+
+
+ for target in selfcal_library:
+ sani_target=sanitize_string(target)
+ for band in selfcal_library[target].keys():
+ if target not in fallback_fields[band]:
+ continue
+ if 'gaintable_final' in selfcal_library[target][band]['vislist'][0]:
+ print('****************Reapplying previous solint solutions*************')
+ for vis in selfcal_library[target][band]['vislist']:
+ print('****************Applying '+str(selfcal_library[target][band][vis]['gaintable_final'])+' to '+target+' '+band+'*************')
+ ## NOTE: should this be selfcal_starting_flags instead of fb_selfcal_starting_flags ???
+ flagmanager(vis=vis,mode='delete',versionname='fb_selfcal_starting_flags_'+sani_target)
+ applycal(vis=vis,\
+ gaintable=selfcal_library[target][band][vis]['gaintable_final'],\
+ interp=selfcal_library[target][band][vis]['applycal_interpolate_final'],\
+ calwt=True,spwmap=selfcal_library[target][band][vis]['spwmap_final'],\
+ applymode=selfcal_library[target][band][vis]['applycal_mode_final'],\
+ field=target,spw=selfcal_library[target][band][vis]['spws'])
+ else:
+ print('****************Removing all calibrations for '+target+' '+band+'**************')
+ for vis in selfcal_library[target][band]['vislist']:
+ flagmanager(vis=vis,mode='delete',versionname='fb_selfcal_starting_flags_'+sani_target)
+ clearcal(vis=vis,field=target,spw=selfcal_library[target][band][vis]['spws'])
+
+ return fallback_fields, calibrators
diff --git a/auto_selfcal/prepare_selfcal.py b/auto_selfcal/prepare_selfcal.py
index 1993cd77..c53b66d2 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'].pop(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'
@@ -534,7 +600,7 @@ def default(self, obj):
return selfcal_library, selfcal_plan, gaincalibrator_dict
-def plan_selfcal_per_solint(selfcal_library, selfcal_plan,optimize_spw_combine=True):
+def plan_selfcal_per_solint(selfcal_library, selfcal_plan, optimize_spw_combine=True, solints=None):
# there are some extra keys in this dictionary that stem from how my thinking was orginally and how it evolved
# the current philosophy is that for each solint it will specify how to apply each gain table, and if it should
# be pre-applied for gaincal solves. Then in gaincal wrapper, the parameters for preapplying all tables and for applying
@@ -559,11 +625,18 @@ def plan_selfcal_per_solint(selfcal_library, selfcal_plan,optimize_spw_combine=T
if selfcal_library[target][band][vis]['baseband'][baseband]['nspws']> maxspws_per_bb:
maxspws_per_bb=selfcal_library[target][band][vis]['baseband'][baseband]['nspws']+0.0
- selfcal_plan[target][band][vis]['solint_settings']={}
- for solint in selfcal_plan[target][band]['solints']:
+ if solints is not None:
+ use_solints = solints
+ else:
+ #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']=[]
@@ -575,27 +648,29 @@ def plan_selfcal_per_solint(selfcal_library, selfcal_plan,optimize_spw_combine=T
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'
selfcal_plan[target][band][vis]['solint_settings'][solint]['final_mode']=''
selfcal_plan[target][band][vis]['solint_settings'][solint]['accepted_gaintable']=''
selfcal_plan[target][band][vis]['solint_settings'][solint]['modes_to_attempt']=[]
- 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])
+
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 solint == 'inf_EB':
+ if "inf_EB" in solint:
selfcal_plan[target][band][vis]['solint_settings'][solint]['modes_to_attempt'].append('combinespwpol')
- selfcal_plan[target][band][vis]['solint_settings'][solint]['preapply_this_gaintable']=True
+ if solint == 'inf_EB':
+ 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']:
+ 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][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
@@ -606,12 +681,15 @@ def plan_selfcal_per_solint(selfcal_library, selfcal_plan,optimize_spw_combine=T
if 'per_bb' not in selfcal_plan[target][band][vis]['solint_settings'][solint]['modes_to_attempt']:
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:
- selfcal_plan[target][band][vis]['solint_settings'][solint]['solmode']='ap'
- else:
- selfcal_plan[target][band][vis]['solint_settings'][solint]['solmode']='p'
- if solint != 'inf_EB' and optimize_spw_combine==False:
- selfcal_plan[target][band][vis]['solint_settings'][solint]['modes_to_attempt']=['combinespw']
+
+ 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'
+
+ if 'inf_EB' not in solint and optimize_spw_combine==False:
+ selfcal_plan[target][band][vis]['solint_settings'][solint]['modes_to_attempt']=['combinespw']
+
for mode in selfcal_plan[target][band][vis]['solint_settings'][solint]['modes_to_attempt']:
gaincal_combine=''
if mode =='combinespw':
@@ -630,11 +708,12 @@ def plan_selfcal_per_solint(selfcal_library, selfcal_plan,optimize_spw_combine=T
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']:
+ 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':
+ 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
selfcal_plan[target][band][vis]['solint_settings'][solint]['filename_append'][mode]=filename_append
@@ -648,7 +727,7 @@ def plan_selfcal_per_solint(selfcal_library, selfcal_plan,optimize_spw_combine=T
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 f0b9165e..fa2bf8f9 100644
--- a/auto_selfcal/run_selfcal.py
+++ b/auto_selfcal/run_selfcal.py
@@ -44,16 +44,26 @@ 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, calibrators[band][0])
+ include_targets, include_scans = triage_calibrators(vislist[0], target, band, calibrators[band][0])
+ print('Co-calibrators: ',include_targets)
+ print('Co-calibrator scans: ',include_scans)
if include_targets == "":
print("No suitable calibrators found, skipping "+target)
selfcal_library['Stop_Reason'] += '; No suitable co-calibrators'
return
-
+ else:
+ for vis in 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'))
+ concatvislist=[]
+ for cal_target in include_targets.split(','):
+ concatvislist.append(vis.replace(sanitize_string(target),sanitize_string(cal_target)))
+ concatvislist.append(vis.replace('.ms','_orig.ms'))
+ print('Concat vislist: ',concatvislist)
+ concat(vis=concatvislist,concatvis=vis)
if selfcal_library['usermodel'] != '':
print('Setting model column to user model')
usermodel_wrapper(selfcal_library,sani_target+'_'+band,
@@ -68,10 +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']):
+ 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'])
+
- print("Solving for solint="+selfcal_plan['solints'][iteration]+' with interval '+selfcal_plan['solint_interval'][iteration])
-
# Set some cocal parameters.
if selfcal_plan['solints'][iteration] in ["inf_EB_fb","inf_fb1"]:
calculate_inf_EB_fb_anyways = True
@@ -81,14 +98,17 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \
preapply_targets_own_inf_EB = False
# If there was not a successful inf_EB solint, then this duplicates inf_fb1 so skip
if "inf_EB" not in selfcal_library[vislist[0]]:
+ iteration += 1
continue
elif not selfcal_library[vislist[0]]["inf_EB"]['Pass']:
+ iteration += 1
continue
elif selfcal_plan['solints'][iteration] == "inf_fb3":
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']:
+ iteration += 1
continue
if 'ap' in selfcal_plan['solints'][iteration]:
@@ -97,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)
+
+ 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]
+ for vis in vislist:
+ 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]
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:
@@ -152,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']
@@ -168,14 +206,13 @@ 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'):
- 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
+ if not checkmask(sani_target+'_'+band+'_'+solint+'_'+str(iteration)+'.image.tt0') and mode !="cocal":
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
@@ -193,8 +230,10 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \
flagmanager(vis=vis, mode = 'restore', versionname = versionname, comment = 'Flag states at start of reduction')
if mode == "cocal":
- flagmanager(vis=vis, mode = 'restore', versionname = 'selfcal_starting_flags', comment = 'Flag states at start of the reduction')
-
+ if os.path.exists(vis+".flagversions/flags.selfcal_starting_flags"):
+ flagmanager(vis=vis, mode = 'restore', versionname = 'selfcal_starting_flags', comment = 'Flag states at start of the reduction')
+ else:
+ flagmanager(vis=vis,mode='save',versionname='selfcal_starting_flags')
if not do_fallback_combinespw:
# We need to redo saving the model now that we have potentially unflagged some data.
if not do_fallback_calonly:
@@ -206,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
@@ -233,38 +280,49 @@ 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,
second_iter_solmode=second_iter_solmode, unflag_fb_to_prev_solint=unflag_fb_to_prev_solint, \
refantmode=refantmode, mode=mode, calibrators=calibrators, gaincalibrator_dict=gaincalibrator_dict,
allow_gain_interpolation=allow_gain_interpolation,spectral_solution_fraction=spectral_solution_fraction,
- do_fallback_calonly=do_fallback_calonly, guess_scan_combine=guess_scan_combine)
+ do_fallback_calonly=do_fallback_calonly, guess_scan_combine=guess_scan_combine, preapply_targets_own_inf_EB=preapply_targets_own_inf_EB)
# 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
@@ -294,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)
@@ -358,7 +419,7 @@ 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) and \
selfcal_plan[fid]['solint_snr_per_field'][solint] > 5)
- 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), \
@@ -436,20 +497,23 @@ 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 (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 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']]):
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 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)
- 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)
-
- 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'
@@ -472,13 +536,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']
@@ -486,8 +550,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]
@@ -496,13 +561,17 @@ 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']:
+ # 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']
@@ -510,8 +579,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
@@ -579,7 +649,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
@@ -592,7 +662,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'
@@ -604,26 +674,38 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \
if reason !='':
reason=reason+'; '
reason=reason+'Beam change beyond '+str(delta_beam_thresh)
- if (post_RMS/RMS > 1.05 and selfcal_plan['solint_snr'][solint] <= 5):
- 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 reason != '':
- reason=reason+'; '
- reason=reason+'NF RMS increase beyond 5%'
+
+ if mode != 'cocal':
+ 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 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%'
+ else:
+ if (post_RMS/RMS > 1.05):
+ if reason != '':
+ reason=reason+'; '
+ reason=reason+'RMS increase beyond 5%'
+ if (post_RMS_NF/RMS_NF > 1.05):
+ if reason != '':
+ reason=reason+'; '
+ reason=reason+'NF RMS increase beyond 5%'
+
if not np.any(field_by_field_success):
if reason != '':
reason=reason+'; '
reason=reason+'All sub-fields failed'
- selfcal_library['Stop_Reason']=reason
for vis in vislist:
+ selfcal_library['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]):
@@ -636,24 +718,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:*************')
@@ -706,25 +788,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 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")
@@ -746,24 +834,34 @@ 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:
print('****************Cocal failed, attempting next inf_fb option*************')
+ iteration += 1
continue
else:
print('****************Aborting further self-calibration attempts for '+target+' '+band+'**************')
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 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 2a7f6da5..69a3e776 100644
--- a/auto_selfcal/selfcal_helpers.py
+++ b/auto_selfcal/selfcal_helpers.py
@@ -516,7 +516,7 @@ def tclean_wrapper(selfcal_library, imagename, band, scales=[0], smallscalebias
uvrange=selfcal_library['uvrange'],
reffreq = reffreq,
threshold=threshold,
- parallel=parallel,
+ parallel=False,
phasecenter=phasecenter,spw=spws_per_vis,wprojplanes=wprojplanes)
elif savemodel=='modelcolumn' and selfcal_library['usermodel'] !='':
@@ -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,42 @@ 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']:
+ 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 +1863,119 @@ 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)]
+ 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
mean_SNR_spw=mean_SNR_spw/total_vis
solint_snr_per_spw[solint][str(spw)]=mean_SNR_spw
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]
+ mean_SNR_bb=mean_SNR_bb/len(vislist)
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=''):
@@ -3058,6 +3077,7 @@ def plot_ants_flagging_colored(filename,vis,gaintable):
def get_flagged_solns_per_ant_from_dict(gc_dict_list,spwlist,vis):
msmd.open(vis)
antids=[]
+ print("len(gc_dict_list)", len(gc_dict_list))
for ant in [idant for idant in gc_dict_list[0]['solvestats']['spw'+str(spwlist[0])].keys() if idant.startswith('ant')]:
antids.append(int(ant.replace('ant','')))
antids.sort()
@@ -3124,6 +3144,8 @@ def plot_ants_flagging_colored_from_dict(filename,selfcal_library,selfcal_plan,s
spwlist_pass=spwlist_bb.copy()
+ print(selfcal_plan.keys())
+ print("plot_ants_flagging_colored_from_dict", solint, final_mode, len(selfcal_plan['solint_settings'][solint]['gaincal_return_dict'][final_mode]))
names, offset_x, offset_y, apriori_flagged, nflagged, nunflagged, ntotal, fracflagged, nflagged_non_apriori, ntotal_non_apriori_flagged, fracflagged_non_apriori=get_flagged_solns_per_ant_from_dict(selfcal_plan['solint_settings'][solint]['gaincal_return_dict'][final_mode],spwlist_pass,vis)
fracflagged=fracflagged_non_apriori
print(fracflagged)
@@ -3628,8 +3650,11 @@ def get_gaincalmode_flagging_stats(selfcal_library,selfcal_plan,vis,gaintable_pr
selfcal_plan[vis]['solint_settings'][solint]['nflags_apriori'][mode],selfcal_plan[vis]['solint_settings'][solint]['nflags'][mode],selfcal_plan[vis]['solint_settings'][solint]['nunflagged'][mode],selfcal_plan[vis]['solint_settings'][solint]['ntotal'][mode],selfcal_plan[vis]['solint_settings'][solint]['fracflagged'][mode],selfcal_plan[vis]['solint_settings'][solint]['nflags_non_apriori'][mode],selfcal_plan[vis]['solint_settings'][solint]['ntotal_non_apriori'][mode],selfcal_plan[vis]['solint_settings'][solint]['fracflagged_non_apriori'][mode]=get_gaintable_flagging_stats(selfcal_plan[vis]['solint_settings'][solint]['gaincal_return_dict'][mode],spwlist_bb)
else:
baseband_scale=1.0
- if solint == 'inf_EB':
+ if 'inf_EB' in solint:
n_solutions=1.0
+ elif 'inf_EB_fb' in selfcal_plan[vis]['solint_settings'].keys():
+ n_antennas=selfcal_plan[vis]['solint_settings']['inf_EB_fb']['ntotal_non_apriori'][coarsest_mode][0]/selfcal_plan[vis]['solint_settings'][solint]['polscale'][mode]
+ n_solutions=(selfcal_plan[vis]['solint_settings'][solint]['nflags_non_apriori'][coarsest_mode][0]+selfcal_plan[vis]['solint_settings'][solint]['nunflagged'][coarsest_mode][0])/n_antennas
elif 'inf_EB' in selfcal_plan[vis]['solint_settings'].keys():
n_antennas=selfcal_plan[vis]['solint_settings']['inf_EB']['ntotal_non_apriori'][coarsest_mode][0]/selfcal_plan[vis]['solint_settings'][solint]['polscale'][mode]
n_solutions=(selfcal_plan[vis]['solint_settings'][solint]['nflags_non_apriori'][coarsest_mode][0]+selfcal_plan[vis]['solint_settings'][solint]['nunflagged'][coarsest_mode][0])/n_antennas
@@ -3765,7 +3790,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'
@@ -3780,7 +3805,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
@@ -4643,11 +4668,23 @@ def unflag_failed_antennas(vis, caltable, gaincal_return, telescope, flagged_fra
-def triage_calibrators(vis, target, potential_calibrators, max_distance=10.0, max_time=600.):
+def triage_calibrators(vis, target, band, potential_calibrators, max_distance=10.0, max_time=600.):
gaincalibrator_dict = {}
+ # account for possible different naming conventions in original visibilities
+ sani_target=sanitize_string(target)
+ orig_vis=''
+ if os.path.exists(vis.replace("_target.selfcal.ms",".ms").replace(sani_target+'_'+band+'_','')):
+ orig_vis=os.path.exists(vis.replace("_target.selfcal.ms",".ms").replace(sani_target+'_'+band+'_'))
+ elif os.path.exists(vis.replace("_targets.selfcal.ms",".ms").replace(sani_target+'_'+band+'_','')):
+ orig_vis=os.path.exists(vis.replace("_targets.selfcal.ms",".ms").replace(sani_target+'_'+band+'_',''))
- if os.path.exists(vis.replace("_target.selfcal.ms",".ms")):
- msmd.open(vis.replace("_target.selfcal.ms",".ms"))
+ # original visibilities, with all sources have a different filename now
+
+ orig_targets_vis=vis.replace(".selfcal.ms",".ms").replace(sani_target+'_'+band+'_','')
+ vis=orig_targets_vis
+
+ if orig_vis !='':
+ msmd.open(orig_vis)
for field in msmd.fieldsforintent("*CALIBRATE_PHASE*"):
scans_for_field = msmd.scansforfield(field)
@@ -4986,9 +5023,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 df02e8ff..c9c06c7d 100644
--- a/auto_selfcal/tests/test_auto_selfcal.py
+++ b/auto_selfcal/tests/test_auto_selfcal.py
@@ -92,6 +92,7 @@ def test_benchmark(tmp_path, dataset):
pytest.param("Band8-7m-2.tar.gz", 'https://nrao-my.sharepoint.com/:u:/g/personal/psheehan_nrao_edu/IQAnFUlx_SR0SLdj_vX8MVUWAbMu9BveHCsb9NB07-OD3bo?e=TNBTXq&download=1', id="Band8-7m-2"),
pytest.param("M82-C-conf-C-band_small.tar.gz", 'https://nrao-my.sharepoint.com/:u:/g/personal/psheehan_nrao_edu/IQCT-CkbAW7YT5LIev5CrDSOAZt5uC_tUIReMwlA14Tu9y4?e=2zF126&download=1', id="M82-C-conf-C-band_small"),
pytest.param("K-band-mini-mosaic.tar.gz", 'https://nrao-my.sharepoint.com/:u:/g/personal/psheehan_nrao_edu/IQBCQc-REEGYQrBwBu2F9uUTAfpCRQq1gYAPO18e-CL_IUk?e=7adSCD&download=1', id='K-band-mini-mosaic'),
+ pytest.param("Band8-7m-cocal.tar.gz", 'https://nrao-my.sharepoint.com/:u:/g/personal/psheehan_nrao_edu/IQAjUGmKxJpoSq25LNp6qnhUAfhcQ-o3RnP7q13r0OgPCfc?e=rX9LaN&download=1', id='Band8-7m-cocal'),
pytest.param("2019.1.00691.S_SB.tar.gz", 'https://nrao-my.sharepoint.com/:u:/g/personal/psheehan_nrao_edu/IQAoeNLrwb1aSrWCfo8ekE6NAXVJ4Qlpps0gdAeY1U9Axn0?e=MMXCaT&download=1', id='2019.1.00691.S_SB'),
pytest.param("VLBA-1-spw.tar.gz", 'https://nrao-my.sharepoint.com/:u:/g/personal/psheehan_nrao_edu/IQATgZY3r_6MS50DqJaFu4QIAVn9jnOq2X-oubxzGtY8-bM?e=fBnvSh&download=1', id='VLBA-1-spw')
]
@@ -121,7 +122,7 @@ def test_on_github(tmp_path, request, zip_file, link):
auto_selfcal(iscalibrator=True, refant='FD,NL,PT', do_delay_cal=True, shorter_amp_solints=True,
targets='J1154+6022', applytargets='J1203+6031', imsize=640, cell='0.0002arcsec')
else:
- auto_selfcal(sort_targets_and_EBs=True, align_EBs=True, weblog=True, usermask=usermask, usermodel=usermodel)
+ auto_selfcal(sort_targets_and_EBs=True, align_EBs=True, weblog=True, allow_cocal=True, usermask=usermask, usermodel=usermodel)
os.system('rm -rf *.ms*') # Delete MS files as space is limited on GitHub.
@@ -130,10 +131,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
@@ -149,13 +167,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={}):
if isinstance(dictionary1, str):
with open(dictionary1, 'rb') as handle:
dictionary1 = pickle.load(handle)
@@ -171,7 +191,7 @@ 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")
else:
@@ -180,20 +200,45 @@ def compare_two_dictionaries(dictionary1, dictionary2, path=[], exclude=[], tole
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
+
+ try:
+ if key not in dictionary1 and key not in key_map and int(key) in dictionary1:
+ key = int(key)
+ except:
+ continue
- if key not in dictionary1 and int(key) in dictionary1:
- key = int(key)
+ 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:
+ print(f"Using alternative key {alt_key} to match with key {key}")
+ else:
+ print(f"No match found in dictionary1, this is a difference")
+ difference_count += 1
+ continue
+ else:
+ alt_key = key
+
- if type(dictionary1[key]) == dict:
- difference_count += compare_two_dictionaries(dictionary1[key], dictionary2[key], path.copy()+[key], exclude=exclude, tolerance=tolerance)
+ 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)
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()
diff --git a/auto_selfcal/weblog_creation.py b/auto_selfcal/weblog_creation.py
index 19e3fff0..61102b7b 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('+'_'+band+'_noise_plot.png)
\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,7 +240,17 @@ 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())
+ 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 | Selfcal stats: |
\n')
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']
for key in quantities:
if key =='Pass':
@@ -270,51 +282,58 @@ 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]):
+ 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 +347,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+=' | \n'
@@ -483,7 +502,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 +512,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 +527,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 +568,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('