diff --git a/auto_selfcal.py b/auto_selfcal.py index e166fec..b5d00f6 100644 --- a/auto_selfcal.py +++ b/auto_selfcal.py @@ -630,7 +630,7 @@ def default(self, obj): selfcal_library[target][band]['per_spw_stats']={} vislist=selfcal_library[target][band]['vislist'].copy() - selfcal_library[target][band]['spw_map'] = get_spw_map(selfcal_library, + selfcal_library[target][band]['spw_map'], selfcal_library[target][band]['reverse_spw_map'] = get_spw_map(selfcal_library, target, band, telescope) #code to work around some VLA data not having the same number of spws due to missing BlBPs @@ -664,6 +664,7 @@ def default(self, obj): for fid in selfcal_library[target][band]['sub-fields']: selfcal_library[target][band][fid]['per_spw_stats']={} selfcal_library[target][band][fid]['spw_map'] = selfcal_library[target][band]['spw_map'] + selfcal_library[target][band][fid]['reverse_spw_map'] = selfcal_library[target][band]['reverse_spw_map'] for vis in selfcal_library[target][band][fid]['vislist']: selfcal_library[target][band][fid][vis]['per_spw_stats'] = {} @@ -779,11 +780,11 @@ def default(self, obj): for solint in solints[band][target]: if solint == 'inf_EB': print('{}: {:0.2f}'.format(solint,solint_snr[target][band][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'])) + print('{}: spw: {}: {:0.2f}'.format(solint,spw,solint_snr_per_spw[target][band][solint][spw])) 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' @@ -873,7 +874,7 @@ def default(self, obj): sourcemodel=usermodel[target] else: sourcemodel='' - run_selfcal(selfcal_library, target, band, solints, solint_snr, solint_snr_per_field, applycal_mode, solmode, band_properties, telescope, n_ants, cellsize[target], imsize[target], \ + run_selfcal(selfcal_library, target, band, solints, solint_snr, solint_snr_per_field, solint_snr_per_spw, applycal_mode, solmode, band_properties, telescope, n_ants, cellsize[target], imsize[target], \ inf_EB_gaintype_dict, inf_EB_gaincal_combine_dict, inf_EB_fallback_mode_dict, gaincal_combine, applycal_interp[target], integration_time, spectral_scan, spws_set,\ gaincal_minsnr=gaincal_minsnr, gaincal_unflag_minsnr=gaincal_unflag_minsnr, minsnr_to_proceed=minsnr_to_proceed, delta_beam_thresh=delta_beam_thresh, do_amp_selfcal=do_amp_selfcal, \ inf_EB_gaincal_combine=inf_EB_gaincal_combine, inf_EB_gaintype=inf_EB_gaintype, unflag_only_lbants=unflag_only_lbants, \ diff --git a/run_selfcal.py b/run_selfcal.py index 751d820..b8d5989 100644 --- a/run_selfcal.py +++ b/run_selfcal.py @@ -12,7 +12,7 @@ except: parallel=False -def run_selfcal(selfcal_library, target, band, solints, solint_snr, solint_snr_per_field, applycal_mode, solmode, band_properties, telescope, n_ants, cellsize, imsize, \ +def run_selfcal(selfcal_library, target, band, solints, solint_snr, solint_snr_per_field, solint_snr_per_spw, applycal_mode, solmode, band_properties, telescope, n_ants, cellsize, imsize, \ inf_EB_gaintype_dict, inf_EB_gaincal_combine_dict, inf_EB_fallback_mode_dict, gaincal_combine, applycal_interp, integration_time, spectral_scan, spws_set, \ gaincal_minsnr=2.0, gaincal_unflag_minsnr=5.0, minsnr_to_proceed=3.0, delta_beam_thresh=0.05, do_amp_selfcal=True, inf_EB_gaincal_combine='scan', inf_EB_gaintype='G', \ unflag_only_lbants=False, unflag_only_lbants_onlyap=False, calonly_max_flagged=0.0, second_iter_solmode="", unflag_fb_to_prev_solint=False, \ @@ -532,7 +532,7 @@ def run_selfcal(selfcal_library, target, band, solints, solint_snr, solint_snr_p solint=solint.replace('_EB','').replace('_ap',''),minsnr=gaincal_minsnr if applymode == "calflag" else max(gaincal_minsnr,gaincal_unflag_minsnr), minblperant=4,combine=test_gaincal_combine, field=include_targets[0],gaintable='',spwmap=[],uvrange=selfcal_library[target][band]['uvrange'], refantmode=refantmode,append=os.path.exists('test_inf_EB_'+gaintype+'.g'))] spwlist=selfcal_library[target][band][vis]['spws'].split(',') - fallback[vis],map_index,spwmap,applycal_spwmap_inf_EB=analyze_inf_EB_flagging(selfcal_library,band,spwlist,sani_target+'_'+vis+'_'+band+'_'+solint+'_'+str(iteration)+'_'+solmode[band][target][iteration]+'.g',vis,target,'test_inf_EB_'+gaincal_gaintype+'.g',spectral_scan,telescope,'test_inf_EB_T.g' if gaincal_gaintype=='G' else None) + fallback[vis],map_index,spwmap,applycal_spwmap_inf_EB=analyze_inf_EB_flagging(selfcal_library,band,spwlist,sani_target+'_'+vis+'_'+band+'_'+solint+'_'+str(iteration)+'_'+solmode[band][target][iteration]+'.g',vis,target,'test_inf_EB_'+gaincal_gaintype+'.g',spectral_scan,telescope, solint_snr_per_spw, minsnr_to_proceed,'test_inf_EB_T.g' if gaincal_gaintype=='G' else None) inf_EB_fallback_mode_dict[target][band][vis]=fallback[vis]+'' print('inf_EB',fallback[vis],applycal_spwmap_inf_EB) diff --git a/selfcal_helpers.py b/selfcal_helpers.py index 3343580..906a641 100644 --- a/selfcal_helpers.py +++ b/selfcal_helpers.py @@ -1184,55 +1184,55 @@ def get_SNR_self_individual(vislist,selfcal_library,n_ant,solints,integration_ti solint_snr = {} solint_snr_per_spw = {} for solint in solints: - #code to work around some VLA data not having the same number of spws due to missing BlBPs - #selects spwlist from the visibilities with the greates number of spws - maxspws=0 - maxspwvis='' - for vis in selfcal_library['vislist']: - if selfcal_library[vis]['n_spws'] >= maxspws: - maxspws=selfcal_library[vis]['n_spws'] - maxspwvis=vis+'' solint_snr[solint]=0.0 solint_snr_per_spw[solint]={} if solint == 'inf_EB': SNR_self_EB=np.zeros(len(selfcal_library['vislist'])) - SNR_self_EB_spw=np.zeros([len(selfcal_library['vislist']),len(selfcal_library[maxspwvis]['spwsarray'])]) - SNR_self_EB_spw_mean=np.zeros([len(selfcal_library[maxspwvis]['spwsarray'])]) SNR_self_EB_spw={} for i in range(len(selfcal_library['vislist'])): SNR_self_EB[i]=SNR/((n_ant)**0.5*(selfcal_library['Total_TOS']/selfcal_library[selfcal_library['vislist'][i]]['TOS'])**0.5) SNR_self_EB_spw[selfcal_library['vislist'][i]]={} - for spw in selfcal_library[selfcal_library['vislist'][i]]['spwsarray']: - if spw in SNR_self_EB_spw[selfcal_library['vislist'][i]].keys(): - SNR_self_EB_spw[selfcal_library['vislist'][i]][str(spw)]=(polscale)**-0.5*SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/selfcal_library[selfcal_library['vislist'][i]]['TOS'])**0.5)*(selfcal_library[selfcal_library['vislist'][i]]['per_spw_stats'][str(spw)]['effective_bandwidth']/selfcal_library[selfcal_library['vislist'][i]]['total_effective_bandwidth'])**0.5 - for spw in selfcal_library[maxspwvis]['spwsarray']: + 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']/selfcal_library[selfcal_library['vislist'][i]]['TOS'])**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 + for spw in selfcal_library['spw_map']: mean_SNR=0.0 + total_vis = 0 for j in range(len(selfcal_library['vislist'])): - if spw in SNR_self_EB_spw[selfcal_library['vislist'][j]].keys(): + if selfcal_library['vislist'][j] in selfcal_library['spw_map'][spw]: mean_SNR+=SNR_self_EB_spw[selfcal_library['vislist'][j]][str(spw)] - mean_SNR=mean_SNR/len(selfcal_library['vislist']) + total_vis += 1 + mean_SNR=mean_SNR/total_vis solint_snr_per_spw[solint][str(spw)]=mean_SNR 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) solint_snr[solint]=selfcal_library['per_scan_SNR'] - for spw in selfcal_library[maxspwvis]['spwsarray']: - 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[maxspwvis]['per_spw_stats'][spw]['effective_bandwidth']/selfcal_library[maxspwvis]['total_effective_bandwidth'])**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']/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 elif solint =='inf' or solint == 'inf_ap': 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[solint]=selfcal_library['per_scan_SNR'] - for spw in selfcal_library[maxspwvis]['spwsarray']: - 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[maxspwvis]['per_spw_stats'][spw]['effective_bandwidth']/selfcal_library[maxspwvis]['total_effective_bandwidth'])**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']/(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 elif solint == 'int': solint_snr[solint]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/integration_time)**0.5) - for spw in selfcal_library[maxspwvis]['spwsarray']: - solint_snr_per_spw[solint][str(spw)]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/integration_time)**0.5)*(selfcal_library[maxspwvis]['per_spw_stats'][spw]['effective_bandwidth']/selfcal_library[maxspwvis]['total_effective_bandwidth'])**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 else: solint_float=float(solint.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[maxspwvis]['spwsarray']: - solint_snr_per_spw[solint][str(spw)]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/solint_float)**0.5)*(selfcal_library[maxspwvis]['per_spw_stats'][spw]['effective_bandwidth']/selfcal_library[maxspwvis]['total_effective_bandwidth'])**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 return solint_snr,solint_snr_per_spw def get_SNR_self_update(all_targets,band,vislist,selfcal_library,n_ant,solint_curr,solint_next,integration_time,solint_snr): @@ -1664,10 +1664,12 @@ def get_spw_map(selfcal_library, target, band, telescope): vislist = [maxspwvis] + vislist spw_map = {} + reverse_spw_map = {} virtual_index = 0 # This code is meant to be generic in order to prepare for cases where multiple EBs might have unique SPWs in them (e.g. inhomogeneous data), # but the criterea for which SPWs match will need to be updated for this to truly generalize. for vis in vislist: + reverse_spw_map[vis] = {} for spw in selfcal_library[target][band][vis]['spwsarray']: found_match = False for s in spw_map: @@ -1705,6 +1707,7 @@ def get_spw_map(selfcal_library, target, band, telescope): if found_match: spw_map[s][vis] = spw + reverse_spw_map[vis][spw] = s break if found_match: @@ -1713,11 +1716,14 @@ def get_spw_map(selfcal_library, target, band, telescope): if not found_match: spw_map[virtual_index] = {} spw_map[virtual_index][vis] = spw + reverse_spw_map[vis][spw] = virtual_index virtual_index += 1 print("spw_map:") print(spw_map) - return spw_map + print("reverse_spw_map:") + print(reverse_spw_map) + return spw_map, reverse_spw_map def largest_prime_factor(n): @@ -3109,7 +3115,7 @@ def get_flagged_solns_per_spw(spwlist,gaintable,extendpol=False): return nflags, nunflagged,fracflagged -def analyze_inf_EB_flagging(selfcal_library,band,spwlist,gaintable,vis,target,spw_combine_test_gaintable,spectral_scan,telescope,spwpol_combine_test_gaintable=None): +def analyze_inf_EB_flagging(selfcal_library,band,spwlist,gaintable,vis,target,spw_combine_test_gaintable,spectral_scan,telescope, solint_snr_per_spw, minsnr_to_proceed,spwpol_combine_test_gaintable=None): if telescope != 'ACA': # if more than two antennas are fully flagged relative to the combinespw results, fallback to combinespw max_flagged_ants_combspw=2.0 @@ -3147,7 +3153,8 @@ def analyze_inf_EB_flagging(selfcal_library,band,spwlist,gaintable,vis,target,sp #if certain spws have more than max_flagged_ants_spwmap flagged solutions that the least flagged spws, set those to spwmap for i in range(len(spwlist)): - if np.min(delta_nflags[i]) > max_flagged_ants_spwmap: + if np.min(delta_nflags[i]) > max_flagged_ants_spwmap or \ + solint_snr_per_spw[target][band]['inf_EB'][str(selfcal_library[target][band]['reverse_spw_map'][vis][int(spwlist[i])])] < minsnr_to_proceed: fallback='spwmap' spwmap[i]=True if total_bws[i] > min_spwmap_bw: