From 507d6b60f47ecc4271a3706f8ae0b0943d86a616 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Mon, 5 Feb 2024 14:34:35 -0500 Subject: [PATCH 01/27] First working attempt at implementing per-EB solints --- gaincal_wrapper.py | 8 +- image_analysis_helpers.py | 7 +- prepare_selfcal.py | 55 +++++++--- run_selfcal.py | 58 +++++++--- selfcal_helpers.py | 218 ++++++++++++++++++++------------------ 5 files changed, 210 insertions(+), 136 deletions(-) diff --git a/gaincal_wrapper.py b/gaincal_wrapper.py index d375c49a..6380ab14 100644 --- a/gaincal_wrapper.py +++ b/gaincal_wrapper.py @@ -260,7 +260,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, ap if mode != 'per_bb': gcdict=gaincal(vis=vis, caltable=gaintable_name, gaintype=selfcal_plan[vis]['solint_settings'][solint]['gaincal_gaintype'], spw=spwselect, refant=selfcal_library[vis]['refant'], calmode=selfcal_plan['solmode'][iteration], solnorm=solnorm if applymode=="calflag" else False, - solint=solint.replace('_EB','').replace('_ap','').replace('scan_','').replace('_fb1','').replace('_fb2','').replace('_fb3',''),\ + solint=selfcal_plan[vis]['solint_settings'][solint]['interval'].replace('_EB','').replace('_ap','').replace('scan_','').replace('_fb1','').replace('_fb2','').replace('_fb3',''),\ minsnr=gaincal_minsnr if applymode == 'calflag' 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',\ @@ -271,7 +271,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, ap spwselect_bb=selfcal_library[vis]['baseband'][baseband]['spwstring'] gcdict=gaincal(vis=vis, caltable=gaintable_name, gaintype=selfcal_plan[vis]['solint_settings'][solint]['gaincal_gaintype'], spw=spwselect_bb, refant=selfcal_library[vis]['refant'], calmode=selfcal_plan['solmode'][iteration], solnorm=solnorm if applymode=="calflag" else False, - solint=solint.replace('_EB','').replace('_ap','').replace('scan_','').replace('_fb1','').replace('_fb2','').replace('_fb3',''),\ + solint=selfcal_plan[vis]['solint_settings'][solint]['interval'].replace('_EB','').replace('_ap','').replace('scan_','').replace('_fb1','').replace('_fb2','').replace('_fb3',''),\ minsnr=gaincal_minsnr if applymode == 'calflag' 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',\ @@ -396,7 +396,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, ap if mode != 'per_bb': gcdict=gaincal(vis=vis, caltable=gaintable_name, gaintype=selfcal_plan[vis]['solint_settings'][solint]['gaincal_gaintype'], spw=spwselect, refant=selfcal_library[vis]['refant'], calmode=selfcal_plan['solmode'][iteration], solnorm=solnorm if applymode=="calflag" else False, - solint=solint.replace('_EB','').replace('_ap','').replace('scan_',''),\ + solint=selfcal_plan[vis]['solint_settings'][solint]['interval'].replace('_EB','').replace('_ap','').replace('scan_',''),\ minsnr=gaincal_minsnr if applymode == 'calflag' 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',\ @@ -407,7 +407,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, ap spwselect_bb=selfcal_library[vis]['baseband'][baseband]['spwstring'] gcdict=gaincal(vis=vis, caltable=gaintable_name, gaintype=selfcal_plan[vis]['solint_settings'][solint]['gaincal_gaintype'], spw=spwselect_bb, refant=selfcal_library[vis]['refant'], calmode=selfcal_plan['solmode'][iteration], solnorm=solnorm if applymode=="calflag" else False, - solint=solint.replace('_EB','').replace('_ap','').replace('scan_',''),\ + solint=selfcal_plan[vis]['solint_settings'][solint]['interval'].replace('_EB','').replace('_ap','').replace('scan_',''),\ minsnr=gaincal_minsnr if applymode == 'calflag' 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',\ diff --git a/image_analysis_helpers.py b/image_analysis_helpers.py index 652bdacf..96bf5593 100644 --- a/image_analysis_helpers.py +++ b/image_analysis_helpers.py @@ -15,7 +15,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/prepare_selfcal.py b/prepare_selfcal.py index eb8f757d..e0411684 100644 --- a/prepare_selfcal.py +++ b/prepare_selfcal.py @@ -446,14 +446,46 @@ 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']=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') + selfcal_plan[target][band] = {} + selfcal_plan[target][band]['solints'] = [] + for vis in selfcal_library[target][band]['vislist']: + selfcal_plan[target][band][vis] = {} + solints,selfcal_plan[target][band][vis]['integration_time'],selfcal_plan[target][band][vis]['gaincal_combine'], \ + selfcal_plan[target][band]['solmode']=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') + + selfcal_plan[target][band][vis]['solint_settings']={} + + subscan_count = 0 + for solint in solints: + """ + 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": + solint_name = "inf_EB" + else: + solint_name = "solint"+str(subscan_count) + if 'ap' in solint: + solint_name += '_ap' + subscan_count += 1 + + if solint_name not in selfcal_plan[target][band]['solints']: + selfcal_plan[target][band]['solints'].append(solint_name) + + selfcal_plan[target][band][vis]['solint_settings'][solint_name]={} + selfcal_plan[target][band][vis]['solint_settings'][solint_name]['interval'] = 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']) - 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']]) ## ## estimate per scan/EB S/N using time on source and median scan times @@ -462,7 +494,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' @@ -493,11 +524,11 @@ 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']={} + #selfcal_plan[target][band][vis]['solint_settings']={} for solint in selfcal_plan[target][band]['solints']: 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']=[] @@ -513,8 +544,8 @@ def plan_selfcal_per_solint(selfcal_library, selfcal_plan,optimize_spw_combine=T selfcal_plan[target][band][vis]['solint_settings'][solint]['accepted_gaintable']='' selfcal_plan[target][band][vis]['solint_settings'][solint]['modes_to_attempt']=[] selfcal_plan[target][band][vis]['solint_settings'][solint]['gaincal_gaintype']='T' - min_SNR_spw=get_min_SNR_spw(selfcal_plan[target][band]['solint_snr_per_spw'][solint]) - min_SNR_bb=get_min_SNR_spw(selfcal_plan[target][band]['solint_snr_per_bb'][solint]) + min_SNR_spw=get_min_SNR_spw(selfcal_plan[target][band][vis]['solint_snr_per_spw'][solint]) + min_SNR_bb=get_min_SNR_spw(selfcal_plan[target][band][vis]['solint_snr_per_bb'][solint]) selfcal_plan[target][band][vis]['solint_settings'][solint]['modes_to_attempt'].append('combinespw') if 'spw' not in selfcal_plan[target][band][vis]['inf_EB_gaincal_combine']: if min_SNR_spw > 2.0: @@ -581,7 +612,7 @@ def set_clean_thresholds(selfcal_library, selfcal_plan, dividing_factor=-99.0, r nsigma_init=np.max([selfcal_library[target][band]['SNR_NF_orig']/dividing_factor_band,5.0]) # count number of amplitude selfcal solints, repeat final clean depth of phase-only for amplitude selfcal - n_ap_solints=sum(1 for solint in selfcal_plan[target][band]['solints'] if 'ap' in solint) + n_ap_solints=sum(1 for solint in selfcal_plan[target][band]['solints'] if 'ap' in selfcal_plan[target][band][selfcal_library[target][band]['vislist'][0]]['solint_settings'][solint]['interval']) if rel_thresh_scaling == 'loge': selfcal_library[target][band]['nsigma'] = np.append(np.exp(np.linspace(np.log(nsigma_init),np.log(3.0),\ diff --git a/run_selfcal.py b/run_selfcal.py index 0e5defa2..b9ee32f0 100644 --- a/run_selfcal.py +++ b/run_selfcal.py @@ -41,8 +41,6 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, telescope, 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]) @@ -56,7 +54,12 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, telescope, n_ants, print('Starting selfcal procedure on: '+target+' '+band) 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]) + for vis in vislist: + print(" "+vis+": "+selfcal_plan['solints'][iteration]) + # Set some cocal parameters. if selfcal_plan['solints'][iteration] in ["inf_EB_fb","inf_fb1"]: @@ -82,11 +85,25 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, telescope, n_ants, break # make sure the last phase-only selfcal table gets pre-applied # solints that use combinespw don't get pre-apply label by default - selfcal_plan[vis]['solint_settings'][selfcal_library['final_phase_solint']]['preapply_this_gaintable']=True + for vis in vislist: + selfcal_plan[vis]['solint_settings'][selfcal_library['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*************') @@ -95,6 +112,11 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, telescope, n_ants, selfcal_library['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] if iteration == 0: print('Starting with solint: '+solint) @@ -479,7 +501,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, telescope, 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 do_fallback_combinespw = False @@ -528,7 +550,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, telescope, n_ants, 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][solint]['Pass']=False selfcal_library[fid][vis][solint]['Fail_Reason']=mosaic_reason[fid] else: @@ -536,7 +558,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, telescope, n_ants, # 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:*************') @@ -569,16 +591,18 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, telescope, n_ants, 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: + 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: + 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. diff --git a/selfcal_helpers.py b/selfcal_helpers.py index 4379d0e4..dab2c23b 100644 --- a/selfcal_helpers.py +++ b/selfcal_helpers.py @@ -988,42 +988,16 @@ 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]['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]['integration_time'], inf_EB_gaincal_combine, - inf_EB_gaintype) - + for vis in selfcal_library[target][band]['vislist']: + 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, selfcal_plan[target][band]['solints'], + 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)) + print(target,band,vis) 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])) + 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'])) @@ -1037,11 +1011,41 @@ 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] = {} + 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, + selfcal_plan[target][band]['solints'], 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 selfcal_plan[target][band]['solints']: + 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,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: @@ -1057,96 +1061,100 @@ def get_SNR_self_individual(vislist,selfcal_library,n_ant,solints,integration_ti #selects spwlist from the visibilities with the greates number of spws maxspws=0 maxspwvis='' - for vis in selfcal_library['vislist']: + #for vis in selfcal_library['vislist']: + for vis in 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]={} solint_snr_per_bb[solint]={} - if solint == 'inf_EB': - SNR_self_EB=np.zeros(len(selfcal_library['vislist'])) + if solint_settings[solint]['interval'] == 'inf_EB': + SNR_self_EB=np.zeros(len(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={} SNR_self_EB_bb={} - 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]]={} - SNR_self_EB_bb[selfcal_library['vislist'][i]]={} - for spw in selfcal_library[selfcal_library['vislist'][i]]['spwsarray']: - 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'][spw]['effective_bandwidth']/selfcal_library[selfcal_library['vislist'][i]]['total_effective_bandwidth'])**0.5 - print(selfcal_library[vis]['baseband']) + for i in range(len(vislist)): + SNR_self_EB[i]=SNR/((n_ant)**0.5*(selfcal_library['Total_TOS']/selfcal_library[vislist[i]]['TOS'])**0.5) + SNR_self_EB_spw[vislist[i]]={} + SNR_self_EB_bb[vislist[i]]={} + for spw in selfcal_library[vislist[i]]['spwsarray']: + SNR_self_EB_spw[vislist[i]][str(spw)]=(polscale)**-0.5*SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/selfcal_library[vislist[i]]['TOS'])**0.5)*(selfcal_library[vislist[i]]['per_spw_stats'][spw]['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[vis]['baseband']: - SNR_self_EB_bb[selfcal_library['vislist'][i]][baseband]=(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]]['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']/selfcal_library[vislist[i]]['TOS'])**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[maxspwvis]['spwsarray']: mean_SNR_spw=0.0 - for j in range(len(selfcal_library['vislist'])): - if str(spw) in SNR_self_EB_spw[selfcal_library['vislist'][j]].keys(): - mean_SNR_spw+=SNR_self_EB_spw[selfcal_library['vislist'][j]][str(spw)] - mean_SNR_spw=mean_SNR_spw/len(selfcal_library['vislist']) + for j in range(len(vislist)): + if str(spw) in SNR_self_EB_spw[vislist[j]].keys(): + mean_SNR_spw+=SNR_self_EB_spw[vislist[j]][str(spw)] + mean_SNR_spw=mean_SNR_spw/len(vislist) print('mean_SNR_spw',mean_SNR_spw,spw) solint_snr_per_spw[solint][str(spw)]=mean_SNR_spw - for baseband in selfcal_library[vis]['baseband']: + for baseband in selfcal_library[maxspwvis]['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]['interval'] =='scan_inf': + selfcal_library['per_scan_SNR']=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/selfcal_library[maxspwvis]['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 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[maxspwvis]['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_per_spw[solint][str(spw)]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/selfcal_library[maxspwvis]['Median_scan_time'])**0.5)*(selfcal_library[maxspwvis]['per_spw_stats'][spw]['effective_bandwidth']/selfcal_library[maxspwvis]['total_effective_bandwidth'])**0.5 + for baseband in selfcal_library[maxspwvis]['baseband']: + solint_snr_per_bb[solint][baseband]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/selfcal_library[maxspwvis]['Median_scan_time'])**0.5)*(selfcal_library[maxspwvis]['baseband'][baseband]['total_effective_bandwidth']/selfcal_library[maxspwvis]['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[maxspwvis]['Median_scan_time']/selfcal_library[maxspwvis]['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 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[maxspwvis]['total_effective_bandwidth'])**0.5 - elif solint == 'int': + solint_snr_per_spw[solint][str(spw)]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/(selfcal_library[maxspwvis]['Median_scan_time']/selfcal_library[maxspwvis]['Median_fields_per_scan']))**0.5)*(selfcal_library[maxspwvis]['per_spw_stats'][spw]['effective_bandwidth']/selfcal_library[maxspwvis]['total_effective_bandwidth'])**0.5 + for baseband in selfcal_library[maxspwvis]['baseband']: + solint_snr_per_bb[solint][baseband]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/(selfcal_library[maxspwvis]['Median_scan_time']/selfcal_library[maxspwvis]['Median_fields_per_scan']))**0.5)*(selfcal_library[maxspwvis]['baseband'][baseband]['total_effective_bandwidth']/selfcal_library[maxspwvis]['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[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 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[maxspwvis]['total_effective_bandwidth'])**0.5 + for baseband in selfcal_library[maxspwvis]['baseband']: + solint_snr_per_bb[solint][baseband]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/integration_time)**0.5)*(selfcal_library[maxspwvis]['baseband'][baseband]['total_effective_bandwidth']/selfcal_library[maxspwvis]['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[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 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[maxspwvis]['total_effective_bandwidth'])**0.5 + for baseband in selfcal_library[maxspwvis]['baseband']: + solint_snr_per_bb[solint][baseband]=SNR/((n_ant-3)**0.5*(selfcal_library['Total_TOS']/solint_float)**0.5)*(selfcal_library[maxspwvis]['baseband'][baseband]['total_effective_bandwidth']/selfcal_library[maxspwvis]['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_ratio=selfcal_library[vis][solint_curr]['SNR_post']/selfcal_library['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[maxspwvis]['spwsarray']: - selfcal_plan['solint_snr_per_spw'][solint_next][str(spw)]=selfcal_plan['solint_snr_per_spw'][solint_next][str(spw)]*SNR_ratio + for spw in selfcal_library[vis]['spwsarray']: + 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 @@ -2538,7 +2546,6 @@ def render_selfcal_solint_summary_table(htmlOut,sclib,target,band,selfcal_plan): line+=''+solint+'\n ' line+='\n' htmlOut.writelines(line) - vis_keys=list(sclib[target][band][vislist[len(vislist)-1]].keys()) quantities=['Pass','intflux_final','intflux_improvement','SNR_final','SNR_Improvement','SNR_NF_final','SNR_NF_Improvement','RMS_final','RMS_Improvement','RMS_NF_final','RMS_NF_Improvement','Beam_Ratio','clean_threshold','Plots'] for key in quantities: if key =='Pass': @@ -2570,51 +2577,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': @@ -2628,7 +2642,7 @@ def render_selfcal_solint_summary_table(htmlOut,sclib,target,band,selfcal_plan): for vis in vislist: line='\n '+vis+': \n' for solint in solint_list: - if solint in vis_keys and sclib[target][band][vis][solint]['Pass'] != 'None' and 'gaintable' in sclib[target][band][vis][solint]: + if solint in sclib[target][band][vis] and sclib[target][band][vis][solint]['Pass'] != 'None' and 'gaintable' in sclib[target][band][vis][solint]: # only evaluate last gaintable not the pre-apply table gaintable=sclib[target][band][vis][solint]['gaintable'][len(sclib[target][band][vis][solint]['gaintable'])-1] line+='antenna positions with flagging plot\n' @@ -2639,7 +2653,7 @@ def render_selfcal_solint_summary_table(htmlOut,sclib,target,band,selfcal_plan): for quantity in ['Nsols_with_preflagged_data','Flagged_Sols_with_preflagged_data','Frac_Flagged_with_preflagged_data','Nsols_without_preflagged_data','Flagged_Sols_without_preflagged_data','Frac_Flagged_without_preflagged_data','SPW_Combine_Mode']: line='\n '+quantity+'\n' for solint in solint_list: - if solint in vis_keys and sclib[target][band][vis][solint]['Pass'] != 'None' and 'gaintable' in sclib[target][band][vis][solint]: + if solint in sclib[target][band][vis] and sclib[target][band][vis][solint]['Pass'] != 'None' and 'gaintable' in sclib[target][band][vis][solint]: # only evaluate last gaintable not the pre-apply table #gaintable=sclib[target][band][vis][solint]['gaintable'][len(sclib[target][band][vis][solint]['gaintable'])-1] #nflagged_sols, nsols=get_sols_flagged_solns(gaintable) @@ -2843,7 +2857,7 @@ def render_per_solint_QA_pages(sclib,selfcal_plan,bands,directory='weblog'): htmlOutSolint.writelines('

Phase vs. Time Plots:

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

MS: '+vis+'

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

') continue ant_list=get_ant_list(vis) From 9a32ff7ce211600a6307657206e7012d70925647 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Tue, 6 Feb 2024 13:30:59 -0500 Subject: [PATCH 02/27] Name the solints p0,p1,p2,...,ap0,ap1,... --- prepare_selfcal.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/prepare_selfcal.py b/prepare_selfcal.py index e0411684..99ba6662 100644 --- a/prepare_selfcal.py +++ b/prepare_selfcal.py @@ -470,10 +470,15 @@ def default(self, obj): if solint == "inf_EB": solint_name = "inf_EB" else: - solint_name = "solint"+str(subscan_count) if 'ap' in solint: - solint_name += '_ap' - subscan_count += 1 + solint_name = 'ap' + 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']: selfcal_plan[target][band]['solints'].append(solint_name) From 4316601dd04ddd99f8f1a91c9cb380cb9ea06b32 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Mon, 29 Jan 2024 14:17:50 -0500 Subject: [PATCH 03/27] Catch cases where spw_baseband == '' in check_spw_widest_in_bb --- selfcal_helpers.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/selfcal_helpers.py b/selfcal_helpers.py index dab2c23b..7f88e9f1 100644 --- a/selfcal_helpers.py +++ b/selfcal_helpers.py @@ -3318,7 +3318,7 @@ def select_best_gaincal_mode(selfcal_library,selfcal_plan,vis,gaintable_prefix,s if np.min(selfcal_plan[vis]['solint_settings'][solint]['delta_nflags']['per_spw'][i]) >= max_flagged_ants_spwmap: fallback='spwmap' spwmap[i]=1.0 - spwmap_widest_window_in_bb=check_spw_widest_in_bb(selfcal_library,vis,spwlist[i]) + spwmap_widest_window_in_bb[i]=check_spw_widest_in_bb(selfcal_library,vis,spwlist[i]) if np.sum(spwmap)/len(spwmap) > 0.5: # if greater than 1/2 of spws need mapping, just assume that we should do combinespw or per_bb fallback='' if np.sum(spwmap_widest_window_in_bb) >= 1.0: # don't do spw mapping within a baseband if the spws to be mapped are the widest in the baseband @@ -3408,6 +3408,8 @@ def find_nearest(array, value): return 1.0 else: return 0.0 + else: + return 1.0 def get_nearest_wide_bw_spw(selfcal_library,vis,spw): mapped_spw=-99 From 5de7f7c4ad48e352dc934f87595ed381e4aa3f49 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Thu, 22 Feb 2024 14:08:36 -0500 Subject: [PATCH 04/27] Fixes to get per-EB solints working when EBs have different #s of solints. --- auto_selfcal.py | 2 +- gaincal_wrapper.py | 8 ++++++-- prepare_selfcal.py | 13 ++++++++++--- run_selfcal.py | 28 +++++++++++++++++----------- selfcal_helpers.py | 39 ++++++++++++++++++++++++--------------- 5 files changed, 58 insertions(+), 32 deletions(-) diff --git a/auto_selfcal.py b/auto_selfcal.py index 697ea706..1166edb9 100644 --- a/auto_selfcal.py +++ b/auto_selfcal.py @@ -524,7 +524,7 @@ def default(self, obj): 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'] line='applycal(vis="'+vis.replace('.selfcal','')+'",gaintable='+str(selfcal_library[target][band][vis]['gaintable_final'])+',interp='+str(selfcal_library[target][band][vis]['applycal_interpolate_final'])+', calwt=False,spwmap='+str(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_orig']+'")\n' applyCalOut.writelines(line) diff --git a/gaincal_wrapper.py b/gaincal_wrapper.py index 6380ab14..017dee4b 100644 --- a/gaincal_wrapper.py +++ b/gaincal_wrapper.py @@ -18,6 +18,8 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, ap ## 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') ## ## Set gaincal parameters depending on which iteration and whether to use combine=spw for inf_EB or not @@ -42,7 +44,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, ap if selfcal_plan['solmode'][iteration]=='p': previous_solint = "inf_EB" else: - previous_solint = selfcal_library['final_phase_solint'] + previous_solint = selfcal_library[vis]['final_phase_solint'] gaincal_spwmap=[] gaincal_preapply_gaintable=[] gaincal_interpolate=[] @@ -77,7 +79,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, ap # 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: + if selfcal_library[vis]['final_solint'] == 'inf_EB' and calculate_inf_EB_fb_anyways: previous_solint = "inf_EB" fallback='' @@ -350,6 +352,8 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, ap 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'] 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']) gaincal_spwmap.append(selfcal_plan[vis]['solint_settings'][selfcal_plan['solints'][j]]['applycal_spwmap']) diff --git a/prepare_selfcal.py b/prepare_selfcal.py index 99ba6662..f8920b82 100644 --- a/prepare_selfcal.py +++ b/prepare_selfcal.py @@ -449,10 +449,11 @@ def default(self, obj): if band in selfcal_library[target]: selfcal_plan[target][band] = {} selfcal_plan[target][band]['solints'] = [] + selfcal_plan[target][band]['solmode'] = [] for vis in selfcal_library[target][band]['vislist']: selfcal_plan[target][band][vis] = {} solints,selfcal_plan[target][band][vis]['integration_time'],selfcal_plan[target][band][vis]['gaincal_combine'], \ - selfcal_plan[target][band]['solmode']=get_solints_simple([vis],scantimesdict[band], + tmp_solmodes=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') @@ -482,6 +483,10 @@ def default(self, obj): if solint_name not in selfcal_plan[target][band]['solints']: selfcal_plan[target][band]['solints'].append(solint_name) + if 'ap' in solint_name: + selfcal_plan[target][band]['solmode'].append('ap') + else: + selfcal_plan[target][band]['solmode'].append('p') selfcal_plan[target][band][vis]['solint_settings'][solint_name]={} selfcal_plan[target][band][vis]['solint_settings'][solint_name]['interval'] = solint @@ -490,7 +495,7 @@ def default(self, obj): print(band,target,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']]) + 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 @@ -531,6 +536,8 @@ def plan_selfcal_per_solint(selfcal_library, selfcal_plan,optimize_spw_combine=T #selfcal_plan[target][band][vis]['solint_settings']={} for solint in selfcal_plan[target][band]['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]={} @@ -617,7 +624,7 @@ def set_clean_thresholds(selfcal_library, selfcal_plan, dividing_factor=-99.0, r nsigma_init=np.max([selfcal_library[target][band]['SNR_NF_orig']/dividing_factor_band,5.0]) # count number of amplitude selfcal solints, repeat final clean depth of phase-only for amplitude selfcal - n_ap_solints=sum(1 for solint in selfcal_plan[target][band]['solints'] if 'ap' in selfcal_plan[target][band][selfcal_library[target][band]['vislist'][0]]['solint_settings'][solint]['interval']) + n_ap_solints=sum(1 for solint in selfcal_plan[target][band]['solints'] if 'ap' in solint) if rel_thresh_scaling == 'loge': selfcal_library[target][band]['nsigma'] = np.append(np.exp(np.linspace(np.log(nsigma_init),np.log(3.0),\ diff --git a/run_selfcal.py b/run_selfcal.py index b9ee32f0..e4425e83 100644 --- a/run_selfcal.py +++ b/run_selfcal.py @@ -58,7 +58,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, telescope, n_ants, print("Solving for solint="+selfcal_plan['solints'][iteration]) for vis in vislist: - print(" "+vis+": "+selfcal_plan['solints'][iteration]) + print(" "+vis+": "+selfcal_plan[vis]['solint_settings'][selfcal_plan['solints'][iteration]]['interval']) # Set some cocal parameters. @@ -86,7 +86,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, telescope, 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": remove_vis = [] @@ -427,8 +427,9 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, telescope, 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] @@ -438,7 +439,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, telescope, n_ants, if field_by_field_success[ind]: selfcal_library[fid]['SC_success']=True selfcal_library[fid]['Stop_Reason']='None' - for vis in selfcal_library[fid]['vislist']: + for vis in selfcal_library[fid]['vislist-to-gaincal']: 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'] @@ -446,8 +447,9 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, telescope, 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 @@ -511,7 +513,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, telescope, 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']: + if not selfcal_library[selfcal_library['vislist-to-gaincal'][0]][solint]['Pass']: reason='' if (post_SNR <= SNR): reason=reason+' S/N decrease' @@ -535,7 +537,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, telescope, n_ants, mosaic_reason = {} new_fields_to_selfcal = [] for fid in selfcal_library['sub-fields-to-selfcal']: - if not selfcal_library[fid][selfcal_library[fid]['vislist'][0]][solint]['Pass']: + if not selfcal_library[fid][selfcal_library[fid]['vislist-to-gaincal'][0]][solint]['Pass']: mosaic_reason[fid]='' if (post_mosaic_SNR[fid] <= mosaic_SNR[fid]): mosaic_reason[fid]=mosaic_reason[fid]+' SNR decrease' @@ -586,12 +588,14 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, telescope, 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') 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']) @@ -599,6 +603,8 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, telescope, n_ants, for fid in selfcal_library['sub-fields-to-selfcal']: 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']) @@ -606,7 +612,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, telescope, n_ants, # 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") diff --git a/selfcal_helpers.py b/selfcal_helpers.py index 7f88e9f1..14078bd0 100644 --- a/selfcal_helpers.py +++ b/selfcal_helpers.py @@ -989,13 +989,15 @@ def get_SNR_self(selfcal_library,selfcal_plan,n_ant,inf_EB_gaincal_combine,inf_E for target in selfcal_library: for band in selfcal_library[target].keys(): 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, selfcal_plan[target][band]['solints'], + 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,vis) - for solint in selfcal_plan[target][band]['solints']: + 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])) ''' @@ -1017,15 +1019,16 @@ def get_SNR_self(selfcal_library,selfcal_plan,n_ant,inf_EB_gaincal_combine,inf_E 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, - selfcal_plan[target][band]['solints'], selfcal_plan[target][band][vis]['solint_settings'], \ + 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 selfcal_plan[target][band]['solints']: + 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])) ''' @@ -2778,7 +2781,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 @@ -2787,7 +2791,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') @@ -2799,7 +2806,7 @@ 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='' for j in range(final_solint_index+index_addition): if selfcal_plan[target][band]['solints'][j] not in keylist: @@ -2813,8 +2820,8 @@ def render_per_solint_QA_pages(sclib,selfcal_plan,bands,directory='weblog'): #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' @@ -2845,17 +2852,17 @@ def render_per_solint_QA_pages(sclib,selfcal_plan,bands,directory='weblog'): htmlOutSolint.writelines('pre-SC-solint image\n') htmlOutSolint.writelines('pre-SC-solint image
\n') - htmlOutSolint.writelines('Post SC SNR: {:0.3f}'.format(sclib[target][band][vislist[0]][selfcal_plan[target][band]['solints'][i]]['SNR_post'])+'
Pre SC SNR: {:0.3f}'.format(sclib[target][band][vislist[0]][selfcal_plan[target][band]['solints'][i]]['SNR_pre'])+'

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

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

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

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

Phase vs. Frequency Plots:

\n') else: htmlOutSolint.writelines('

Phase vs. Time Plots:

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

MS: '+vis+'

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

') @@ -4106,10 +4113,12 @@ def get_min_SNR_spw(snr_per_spw): return minsnr def remove_modes(selfcal_plan,vis,start_index): + 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 - preferred_mode=selfcal_plan[vis]['solint_settings'][selfcal_plan['solints'][j]]['final_mode'] 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']: selfcal_plan[vis]['solint_settings'][selfcal_plan['solints'][j]]['modes_to_attempt'].remove('per_spw') From 7f477232ddcb4965de3f454e3f2d9d9e9e169f54 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Tue, 1 Apr 2025 06:58:11 -0500 Subject: [PATCH 05/27] Show solution interval per-EB in the weblog. --- auto_selfcal/weblog_creation.py | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/auto_selfcal/weblog_creation.py b/auto_selfcal/weblog_creation.py index 0fa4439c..9d112468 100644 --- a/auto_selfcal/weblog_creation.py +++ b/auto_selfcal/weblog_creation.py @@ -221,6 +221,17 @@ def render_selfcal_solint_summary_table(htmlOut,sclib,target,band,selfcal_plan): line+=''+solint+'\n ' line+='\n' htmlOut.writelines(line) + 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': @@ -492,13 +503,20 @@ def render_per_solint_QA_pages(sclib,selfcal_plan,bands,directory='weblog'): htmlOutSolint.writelines('

Targets:

\n') #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') From b22339ab7789459a8ed00b98d0132a86712da292 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Tue, 8 Jul 2025 18:34:55 +0000 Subject: [PATCH 06/27] Add option to use uniform solution intervals across all EBs --- auto_selfcal/__main__.py | 1 + auto_selfcal/auto_selfcal.py | 3 ++- auto_selfcal/prepare_selfcal.py | 21 +++++++++++++++++---- 3 files changed, 20 insertions(+), 5 deletions(-) diff --git a/auto_selfcal/__main__.py b/auto_selfcal/__main__.py index 7f6b8237..87c3f7d7 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 52415825..2f62b2b6 100644 --- a/auto_selfcal/auto_selfcal.py +++ b/auto_selfcal/auto_selfcal.py @@ -46,6 +46,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, @@ -134,7 +135,7 @@ def auto_selfcal( selfcal_library, selfcal_plan, gaincalibrator_dict = prepare_selfcal(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,debug=debug) + uniform_solints=uniform_solints, usermask=usermask, usermodel=usermodel,debug=debug) with open('selfcal_library.pickle', 'wb') as handle: diff --git a/auto_selfcal/prepare_selfcal.py b/auto_selfcal/prepare_selfcal.py index 337e73c9..bd0211af 100644 --- a/auto_selfcal/prepare_selfcal.py +++ b/auto_selfcal/prepare_selfcal.py @@ -9,6 +9,7 @@ def prepare_selfcal(vislist, inf_EB_gaintype='G', apply_cal_mode_default='calflag', do_amp_selfcal=True, + uniform_solints=False, usermask={}, usermodel={}, debug=False): @@ -492,13 +493,25 @@ def default(self, obj): selfcal_plan[target][band] = {} selfcal_plan[target][band]['solints'] = [] selfcal_plan[target][band]['solmode'] = [] - for vis in selfcal_library[target][band]['vislist']: - selfcal_plan[target][band][vis] = {} - solints,selfcal_plan[target][band][vis]['integration_time'],selfcal_plan[target][band][vis]['gaincal_combine'], \ - tmp_solmodes=get_solints_simple([vis],scantimesdict[band], + + if uniform_solints: + solints,tmp_integration_time,tmp_gaincal_combine, \ + tmp_solmodes=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') + 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=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') + 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 From 4d9b729388019e39b2c1627324573677ba06f0e9 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Tue, 8 Jul 2025 18:56:44 +0000 Subject: [PATCH 07/27] Use uniform_solints=True for tests for now to match what was previously happening so that tests pass (if changes haven't messed that up). --- auto_selfcal/tests/test_auto_selfcal.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/auto_selfcal/tests/test_auto_selfcal.py b/auto_selfcal/tests/test_auto_selfcal.py index 04ab1f8e..3b292ec7 100644 --- a/auto_selfcal/tests/test_auto_selfcal.py +++ b/auto_selfcal/tests/test_auto_selfcal.py @@ -50,7 +50,7 @@ def test_benchmark(tmp_path, dataset): ex.update_parameters(partition="batch2", nodes=1, ntasks_per_node=8, cpus_per_task=1, use_srun=False, time=10080, \ mem="128gb", job_name=dataset) - job = ex.submit(auto_selfcal, sort_targets_and_EBs=True, weblog=True, parallel=True) + job = ex.submit(auto_selfcal, sort_targets_and_EBs=True, uniform_solints=True, weblog=True, parallel=True) job.wait() assert job.state in ['DONE','COMPLETED'] @@ -83,7 +83,7 @@ def test_on_github(tmp_path, request, zip_file, link): os.system(f'tar xf {zip_file}') os.system(f'rm -rf {zip_file}') - auto_selfcal(sort_targets_and_EBs=True, weblog=True) + auto_selfcal(sort_targets_and_EBs=True, uniform_solints=True, weblog=True) os.system('rm -rf *.ms*') # Delete MS files as space is limited on GitHub. From 1562607988b995a3d72fa1052dc1f87ee2eeb9da Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Thu, 25 Sep 2025 16:42:53 +0000 Subject: [PATCH 08/27] Add code to enable testing across the major changes brought by per_EB_solints. --- .github/workflows/run_E2E_test.yml | 2 +- auto_selfcal/tests/test_auto_selfcal.py | 44 ++++++++++++++++++++----- 2 files changed, 36 insertions(+), 10 deletions(-) diff --git a/.github/workflows/run_E2E_test.yml b/.github/workflows/run_E2E_test.yml index c02f4852..fd1f0d40 100644 --- a/.github/workflows/run_E2E_test.yml +++ b/.github/workflows/run_E2E_test.yml @@ -5,7 +5,7 @@ name: E2E Tests on: push: - branches: [ "make_pip_installable" ] + branches: [ "make_pip_installable", "make_pip_installable+per_EB_solints+testing" ] pull_request: branches: [ "make_pip_installable" ] diff --git a/auto_selfcal/tests/test_auto_selfcal.py b/auto_selfcal/tests/test_auto_selfcal.py index 3b292ec7..ce14a783 100644 --- a/auto_selfcal/tests/test_auto_selfcal.py +++ b/auto_selfcal/tests/test_auto_selfcal.py @@ -92,7 +92,20 @@ 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) + 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] = [] + + solint_map[solint].append(selfcal_plan[target][band][vis]['solint_settings'][solint]['interval']) + + difference_count = compare_two_dictionaries(selfcal_library1, selfcal_library2, tolerance=0.001, key_map=solint_map) assert difference_count == 0 @@ -114,7 +127,7 @@ def compare_values(list1, list2, tol=1e-3): 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={}): difference_count = 0 all_keys = np.unique(list(dictionary1.keys()) + list(dictionary2.keys())) @@ -123,7 +136,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: @@ -132,20 +145,33 @@ 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]): + continue - if key not in dictionary1 and int(key) in dictionary1: - key = int(key) + try: + if key not in dictionary1 and int(key) in dictionary1: + key = int(key) + except: + continue + + if key in dictionary2 and not key in dictionary1 and key in key_map: + for alt_key in key_map: + if alt_key in dictionary1: + break + 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() From 15edeea62d2e86cd45ca2d64bd08d7f3e430492e Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Tue, 30 Sep 2025 19:20:31 +0000 Subject: [PATCH 09/27] Make sure the SNR_self_EB array is only as long as the vislist passed to the function. --- auto_selfcal/selfcal_helpers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/auto_selfcal/selfcal_helpers.py b/auto_selfcal/selfcal_helpers.py index cf5e2ad7..f049b925 100644 --- a/auto_selfcal/selfcal_helpers.py +++ b/auto_selfcal/selfcal_helpers.py @@ -1283,7 +1283,7 @@ def get_SNR_self_individual(vislist,selfcal_library,n_ant,solints,solint_setting solint_snr_per_spw[solint]={} solint_snr_per_bb[solint]={} if solint_settings[solint]['interval'] == 'inf_EB': - 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(vislist)): From 3b50cb0807ce9b02e8f4c99124adaf80379f39da Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Tue, 30 Sep 2025 19:23:20 +0000 Subject: [PATCH 10/27] Add some keywords to exclude from comparison due to expected changes in this branch. --- auto_selfcal/tests/test_auto_selfcal.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/auto_selfcal/tests/test_auto_selfcal.py b/auto_selfcal/tests/test_auto_selfcal.py index ce14a783..b8d91a21 100644 --- a/auto_selfcal/tests/test_auto_selfcal.py +++ b/auto_selfcal/tests/test_auto_selfcal.py @@ -105,7 +105,8 @@ def test_on_github(tmp_path, request, zip_file, link): solint_map[solint].append(selfcal_plan[target][band][vis]['solint_settings'][solint]['interval']) - difference_count = compare_two_dictionaries(selfcal_library1, selfcal_library2, tolerance=0.001, key_map=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"]) assert difference_count == 0 From d23ce03c1917a7af01c1db939af9d5292fd7a358 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Wed, 8 Apr 2026 18:35:00 +0000 Subject: [PATCH 11/27] Ensure inf_EB and *_ap solints are mapped correctly now that we record the interval instead of the actual stage name in the solint settings. Also enable proper comparison of None types. --- auto_selfcal/tests/test_auto_selfcal.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/auto_selfcal/tests/test_auto_selfcal.py b/auto_selfcal/tests/test_auto_selfcal.py index f880cf0c..179a4d10 100644 --- a/auto_selfcal/tests/test_auto_selfcal.py +++ b/auto_selfcal/tests/test_auto_selfcal.py @@ -116,7 +116,14 @@ def test_on_github(tmp_path, request, zip_file, link): if solint not in solint_map: solint_map[solint] = [] - solint_map[solint].append(selfcal_plan[target][band][vis]['solint_settings'][solint]['interval']) + mapped_solint = selfcal_plan[target][band][vis]['solint_settings'][solint]['interval'] + if solint == 'inf_EB': + mapped_solint += '_EB' + elif 'ap' in solint: + mapped_solint += '_ap' + + 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"]) @@ -135,6 +142,8 @@ 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 From c28ed0902d94da9bed51c3fc55dce98c27817c71 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Fri, 24 Apr 2026 20:13:31 +0000 Subject: [PATCH 12/27] Properly handle the key_map so that keys whose names changed are looked up properly. Also add some informative messaging, and add telescope to the list of excluded keys because it is new. --- auto_selfcal/tests/test_auto_selfcal.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/auto_selfcal/tests/test_auto_selfcal.py b/auto_selfcal/tests/test_auto_selfcal.py index abba2447..746ced34 100644 --- a/auto_selfcal/tests/test_auto_selfcal.py +++ b/auto_selfcal/tests/test_auto_selfcal.py @@ -130,7 +130,7 @@ def test_on_github(tmp_path, request, zip_file, link): 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"]) + exclude=["final_phase_solint", "final_solint", "gaintable_final", "per_EB_SNR", "vislist-to-gaincal", "telescope"]) assert difference_count == 0 @@ -180,18 +180,30 @@ def compare_two_dictionaries(dictionary1, dictionary2, path=[], exclude=[], tole 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 int(key) in dictionary1: + if key not in dictionary1 and key not in key_map and int(key) in dictionary1: key = int(key) except: continue if key in dictionary2 and not key in dictionary1 and key in key_map: - for alt_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 From 38f5711a8385376c23f620115c2171729646195e Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Thu, 7 May 2026 19:42:11 +0000 Subject: [PATCH 13/27] Track the sub-fields to selfcal/gaincal per-EB. Also store the original solint name in the solint settings to keep track of special solints like scan_inf and 300s_ap --- auto_selfcal/gaincal_wrapper.py | 18 +++--- auto_selfcal/mosaic_helpers.py | 69 +++++++++++---------- auto_selfcal/prepare_selfcal.py | 9 ++- auto_selfcal/run_selfcal.py | 81 ++++++++++++++++++------- auto_selfcal/selfcal_helpers.py | 4 +- auto_selfcal/tests/test_auto_selfcal.py | 7 ++- 6 files changed, 117 insertions(+), 71 deletions(-) diff --git a/auto_selfcal/gaincal_wrapper.py b/auto_selfcal/gaincal_wrapper.py index 3a9da755..7a3b3bfd 100644 --- a/auto_selfcal/gaincal_wrapper.py +++ b/auto_selfcal/gaincal_wrapper.py @@ -109,7 +109,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, ap include_targets = selfcal_library['bands_for_targets'][vis]['field_str'] include_scans = "" - if solint == "scan_inf": + if selfcal_plan[vis]['solint_settings'][solint]['sub-name'] == "scan_inf": if len(gaincalibrator_dict[vis]) > 0: print("Determining scan_inf from calibrator scans in full MS") scans = [] @@ -285,10 +285,10 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, ap for incl_scan in include_scans: scan_targets = [] for fid in [selfcal_library['sub-fields-fid_map'][vis][fid] for fid in \ - np.intersect1d(selfcal_library['sub-fields-to-gaincal'],list(selfcal_library['sub-fields-fid_map'][vis].keys()))] if incl_scan == '' else \ + np.intersect1d(selfcal_library[vis]['sub-fields-to-gaincal'],list(selfcal_library['sub-fields-fid_map'][vis].keys()))] if incl_scan == '' else \ np.intersect1d(msmd.fieldsforscans(np.array(incl_scan.split(",")).astype(int)), \ [selfcal_library['sub-fields-fid_map'][vis][fid] for fid in \ - numpy.intersect1d(selfcal_library['sub-fields-to-gaincal'],list(selfcal_library['sub-fields-fid_map'][vis].keys()))]): + numpy.intersect1d(selfcal_library[vis]['sub-fields-to-gaincal'],list(selfcal_library['sub-fields-fid_map'][vis].keys()))]): # Note: because of the msmd above getting actual fids from the MS, we just need to append fid below. scan_targets.append(fid) @@ -412,7 +412,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, ap 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]+'' @@ -430,7 +430,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, ap 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=[] @@ -537,7 +537,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, ap 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]+'' @@ -612,7 +612,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, ap 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. @@ -656,7 +656,7 @@ def gaincal_wrapper(selfcal_library, selfcal_plan, target, band, vis, solint, ap 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() @@ -714,7 +714,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/mosaic_helpers.py b/auto_selfcal/mosaic_helpers.py index 93543f2a..3301bfde 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]['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][solint]['Pass'] = 'None' + if 'Fail_Reason' in selfcal_library[fid][vis][solint]: + selfcal_library[fid][vis][solint]['Fail_Reason'] += '; ' + else: + selfcal_library[fid][vis][solint]['Fail_Reason'] = '' + selfcal_library[fid][vis][solint]['Fail_Reason'] += 'No viable fields' + return new_fields_to_selfcal ## NEXT TO DO: check % of flagged solutions - DONE, see above ## After that enable option for interpolation through inf - DONE tb.open(sani_target+'_'+vis+'_'+band+'_'+solint+'_'+str(iteration)+'_'+solmode[iteration]+'_'+ @@ -100,7 +99,7 @@ def evaluate_subfields_after_gaincal(selfcal_library, target, band, solint, iter scans = tb.getcol("SCAN_NUMBER") for fid in np.intersect1d(new_fields_to_selfcal,list(selfcal_library['sub-fields-fid_map'][vis].keys())): - if solint == "scan_inf": + if selfcal_plan[vis]['solint_settings'][solint]['sub-name'] == "scan_inf": msmd.open(vis) cals_for_scan = [] total_cals_for_scan = [] @@ -125,29 +124,29 @@ def evaluate_subfields_after_gaincal(selfcal_library, target, band, solint, iter new_fields_to_selfcal.remove(fid) if fid not in new_fields_to_selfcal: - # We need to update all the EBs, not just the one that failed. - for v in selfcal_library[fid]['vislist']: - selfcal_library[fid][v][solint]['Pass'] = 'None' - if allow_gain_interpolation: - selfcal_library[fid][v][solint]['Fail_Reason'] = 'Interpolation beyond inf' - else: - selfcal_library[fid][v][solint]['Fail_Reason'] = 'Bad gaincal solutions' + selfcal_library[fid][vis][solint]['Pass'] = 'None' + if allow_gain_interpolation: + selfcal_library[fid][vis][solint]['Fail_Reason'] = 'Interpolation beyond inf' + else: + selfcal_library[fid][vis][solint]['Fail_Reason'] = 'Bad gaincal solutions' tb.close() elif selfcal_library['obstype'] == 'mosaic' and solint == "inf_EB": ## If an EB had no fields to gaincal on, remove all fields in that EB from being selfcal'd as there is no calibration available ## in this EB. - for vis in selfcal_library['vislist']: - if np.intersect1d(selfcal_library['sub-fields-to-gaincal'],\ + #for vis in selfcal_library['vislist']: + if True: + if np.intersect1d(selfcal_library[vis]['sub-fields-to-gaincal'],\ list(selfcal_library['sub-fields-fid_map'][vis].keys())).size == 0: for fid in np.intersect1d(new_fields_to_selfcal,list(selfcal_library['sub-fields-fid_map'][vis].keys())): new_fields_to_selfcal.remove(fid) selfcal_library[fid]['Stop_Reason'] = 'No viable calibrator fields for inf_EB in at least 1 EB' - for v in selfcal_library[fid]['vislist']: - selfcal_library[fid][v][solint]['Pass'] = 'None' - selfcal_library[fid][v][solint]['Fail_Reason'] = 'No viable inf_EB fields' + selfcal_library[fid][vis][solint]['Pass'] = 'None' + selfcal_library[fid][vis][solint]['Fail_Reason'] = 'No viable inf_EB fields' + + print("new_fields_to_selfcal", new_fields_to_selfcal) return new_fields_to_selfcal diff --git a/auto_selfcal/prepare_selfcal.py b/auto_selfcal/prepare_selfcal.py index b294397a..1cb94dad 100644 --- a/auto_selfcal/prepare_selfcal.py +++ b/auto_selfcal/prepare_selfcal.py @@ -152,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. @@ -161,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] = {} @@ -568,6 +572,7 @@ def default(self, obj): 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']) @@ -690,11 +695,11 @@ 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']: 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': 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 diff --git a/auto_selfcal/run_selfcal.py b/auto_selfcal/run_selfcal.py index 67625997..0a3a00dc 100644 --- a/auto_selfcal/run_selfcal.py +++ b/auto_selfcal/run_selfcal.py @@ -66,6 +66,9 @@ 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']] @@ -172,7 +175,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'] @@ -195,7 +198,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ for vis in vislist: 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][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 @@ -226,13 +229,22 @@ 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['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] @@ -253,13 +265,17 @@ 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 @@ -275,16 +291,23 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ # With gaincal done and bad fields removed from gain tables if necessary, check whether any fields should no longer be # selfcal'd because they have too much interpolation. if selfcal_library['obstype'] == 'mosaic': - selfcal_library['sub-fields-to-selfcal'] = evaluate_subfields_after_gaincal(selfcal_library, target, band, - solint, iteration, selfcal_plan['solmode'], allow_gain_interpolation=allow_gain_interpolation) + selfcal_library['sub-fields-to-selfcal'] = [] + for vis in vislist: + selfcal_library[vis]['sub-fields-to-selfcal'] = evaluate_subfields_after_gaincal(vis, selfcal_library, selfcal_plan, target, band, + solint, iteration, selfcal_plan['solmode'], allow_gain_interpolation=allow_gain_interpolation) + selfcal_library['sub-fields-to-selfcal'] = np.union1d(selfcal_library['sub-fields-to-selfcal'], selfcal_library[vis]['sub-fields-to-selfcal']).astype(int).tolist() + + print('selfcal_library["sub-fields-to-selfcal"] = ', selfcal_library['sub-fields-to-selfcal']) + for vis in selfcal_library['vislist']: + print('selfcal_library[vis]["sub-fields-to-selfcal"] = ', selfcal_library[vis]['sub-fields-to-selfcal']) ## ## Apply gain solutions per MS, target, solint, and band ## for vis in vislist: applycal_wrapper(vis, target, band, solint, selfcal_library, - current=lambda f: f in selfcal_library['sub-fields-to-selfcal'], - final=lambda f: f not in selfcal_library['sub-fields-to-selfcal'] and selfcal_library[f]['SC_success'], + current=lambda f: f in selfcal_library[vis]['sub-fields-to-selfcal'], + final=lambda f: f not in selfcal_library[vis]['sub-fields-to-selfcal'] and selfcal_library[f]['SC_success'], restore_flags='fb_selfcal_starting_flags_'+sani_target+'_'+band if mode == "cocal" else None) ## Create post self-cal image using the model as a startmodel to evaluate how much selfcal helped @@ -314,6 +337,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) @@ -524,6 +550,10 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ selfcal_library[fid]['inf_EB_SNR_decrease']=False 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]['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'] @@ -774,6 +804,12 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ if mode == "selfcal" and iteration > 1 and selfcal_plan['solmode'][iteration] !='ap' and do_amp_selfcal: # if a solution interval shorter than inf for phase-only SC has passed, attempt amplitude selfcal iteration=selfcal_plan['solmode'].index('ap') selfcal_library['sub-fields-to-selfcal'] = selfcal_library['sub-fields'] + for vis in vislist: + selfcal_library[vis]['sub-fields-to-selfcal'] = selfcal_library['sub-fields'].copy() + for fid in selfcal_library[vis]['sub-fields-to-selfcal']: + if not fid in selfcal_library['sub-fields-fid_map'][vis]: + # If a sub-field is not in an EB, it shouldn't be considered for selfcal for that EB + selfcal_library[vis]['sub-fields-to-selfcal'].pop(fid) print('****************Selfcal halted for phase, attempting amplitude*************') continue elif mode == "cocal" and "inf_fb" in solint: @@ -784,14 +820,17 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ break # breakout of loops of successive solints since solutions are getting worse # Finally, update the list of fields to be self-calibrated now that we don't need to know the list at the beginning of this solint. - new_fields_to_selfcal = [] coarsest_solint=selfcal_plan['solints'][0] - for fid in selfcal_library['sub-fields']: - if mode == "cocal": - if ("inf_EB" in selfcal_library[fid]['vislist'][0] and selfcal_library[fid][selfcal_library[fid]['vislist'][0]][coarsest_solint]["Pass"]) or selfcal_library[fid][selfcal_library[fid]['vislist'][0]][coarsest_solint+"_fb"]["Pass"]: - new_fields_to_selfcal.append(fid) - else: - if selfcal_library[fid][selfcal_library[fid]['vislist'][0]][coarsest_solint]["Pass"]: - new_fields_to_selfcal.append(fid) + selfcal_library['sub-fields-to-selfcal'] = [] + for vis in vislist: + new_fields_to_selfcal = [] + for fid in selfcal_library['sub-fields']: + if 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() \ No newline at end of file diff --git a/auto_selfcal/selfcal_helpers.py b/auto_selfcal/selfcal_helpers.py index 6fa9b89c..210a2c18 100644 --- a/auto_selfcal/selfcal_helpers.py +++ b/auto_selfcal/selfcal_helpers.py @@ -1906,7 +1906,7 @@ def get_SNR_self_individual(vislist,selfcal_library,n_ant,solints,solint_setting 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_settings[solint]['interval'] =='scan_inf': + 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']: @@ -3799,7 +3799,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 diff --git a/auto_selfcal/tests/test_auto_selfcal.py b/auto_selfcal/tests/test_auto_selfcal.py index 746ced34..d4c9eca9 100644 --- a/auto_selfcal/tests/test_auto_selfcal.py +++ b/auto_selfcal/tests/test_auto_selfcal.py @@ -120,7 +120,7 @@ def test_on_github(tmp_path, request, zip_file, link): if solint not in solint_map: solint_map[solint] = [] - mapped_solint = selfcal_plan[target][band][vis]['solint_settings'][solint]['interval'] + mapped_solint = selfcal_plan[target][band][vis]['solint_settings'][solint]['sub-name'] if solint == 'inf_EB': mapped_solint += '_EB' elif 'ap' in solint: @@ -130,7 +130,10 @@ def test_on_github(tmp_path, request, zip_file, link): 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"]) + 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"]) assert difference_count == 0 From 4392918a2b422cf0e7068480ed8abeafabf302c8 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Mon, 11 May 2026 18:30:26 +0000 Subject: [PATCH 14/27] No need to add _ap or _EB for solint matching now that we are tracking the sub-name. --- auto_selfcal/tests/test_auto_selfcal.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/auto_selfcal/tests/test_auto_selfcal.py b/auto_selfcal/tests/test_auto_selfcal.py index d4c9eca9..2fba8761 100644 --- a/auto_selfcal/tests/test_auto_selfcal.py +++ b/auto_selfcal/tests/test_auto_selfcal.py @@ -121,10 +121,6 @@ def test_on_github(tmp_path, request, zip_file, link): solint_map[solint] = [] mapped_solint = selfcal_plan[target][band][vis]['solint_settings'][solint]['sub-name'] - if solint == 'inf_EB': - mapped_solint += '_EB' - elif 'ap' in solint: - mapped_solint += '_ap' solint_map[solint].append(mapped_solint) print(solint_map) From e733a4aed5a7516e9bab9f4dedcb23d3b1f8cc64 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Tue, 12 May 2026 09:22:23 -0400 Subject: [PATCH 15/27] Make sure the ap solints actually use solmode ap --- auto_selfcal/prepare_selfcal.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/auto_selfcal/prepare_selfcal.py b/auto_selfcal/prepare_selfcal.py index 1cb94dad..ff1e7e64 100644 --- a/auto_selfcal/prepare_selfcal.py +++ b/auto_selfcal/prepare_selfcal.py @@ -671,7 +671,7 @@ 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: + 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' From 5fe98a6a11f5fd75ab8caeaa05504f6c8edf8239 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Tue, 12 May 2026 10:03:40 -0400 Subject: [PATCH 16/27] sub-fields-to-selfcal was not being properly updated per-vis after evaluating which fields should be gaincal'd --- auto_selfcal/run_selfcal.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/auto_selfcal/run_selfcal.py b/auto_selfcal/run_selfcal.py index 0a3a00dc..4ce4eb33 100644 --- a/auto_selfcal/run_selfcal.py +++ b/auto_selfcal/run_selfcal.py @@ -239,7 +239,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ if solint != 'inf_EB' and not allow_gain_interpolation: selfcal_library['sub-fields-to-selfcal'] = selfcal_library['sub-fields-to-gaincal'] for vis in vislist: - selfcal_library['sub-fields-to-selfcal'] = selfcal_library['sub-fields-to-gaincal'] + selfcal_library[vis]['sub-fields-to-selfcal'] = selfcal_library['sub-fields-to-gaincal'] print('Fields to gaincal: ') for vis in vislist: @@ -833,4 +833,4 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ new_fields_to_selfcal.append(fid) 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() \ No newline at end of file + selfcal_library['sub-fields-to-selfcal'] = np.union1d(selfcal_library['sub-fields-to-selfcal'], selfcal_library[vis]['sub-fields-to-selfcal']).tolist() From 7da3c3f5e5213640820555357d37f382595e92eb Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Tue, 12 May 2026 17:46:58 +0000 Subject: [PATCH 17/27] Track Stop_Reason per EB as well. --- auto_selfcal/mosaic_helpers.py | 6 ++--- auto_selfcal/run_selfcal.py | 29 ++++++++++++------------- auto_selfcal/tests/test_auto_selfcal.py | 2 +- auto_selfcal/weblog_creation.py | 16 ++++++++------ 4 files changed, 27 insertions(+), 26 deletions(-) diff --git a/auto_selfcal/mosaic_helpers.py b/auto_selfcal/mosaic_helpers.py index 3301bfde..65a470a4 100644 --- a/auto_selfcal/mosaic_helpers.py +++ b/auto_selfcal/mosaic_helpers.py @@ -55,7 +55,7 @@ def evaluate_subfields_to_gaincal(vis, selfcal_library, target, band, solint, it if fid not in new_fields_to_selfcal and solint != "inf_EB" and not allow_gain_interpolation: #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]['Stop_Reason'] = skip_reason selfcal_library[fid][vis][solint]['Pass'] = "None" selfcal_library[fid][vis][solint]['Fail_Reason'] = skip_reason @@ -83,7 +83,7 @@ def evaluate_subfields_after_gaincal(vis, selfcal_library, selfcal_plan, target, 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' + 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'] += '; ' @@ -142,7 +142,7 @@ def evaluate_subfields_after_gaincal(vis, selfcal_library, selfcal_plan, target, 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' + 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' diff --git a/auto_selfcal/run_selfcal.py b/auto_selfcal/run_selfcal.py index 4ce4eb33..f69bbbed 100644 --- a/auto_selfcal/run_selfcal.py +++ b/auto_selfcal/run_selfcal.py @@ -126,9 +126,10 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ 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 @@ -191,14 +192,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'): 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[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 @@ -246,17 +246,16 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ 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 @@ -518,13 +517,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 '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 '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'] @@ -543,7 +542,6 @@ 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]: @@ -554,6 +552,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ # 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'] @@ -668,8 +667,8 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ if reason != '': reason=reason+'; ' reason=reason+'All sub-fields failed' - selfcal_library['Stop_Reason']=reason for vis in vislist: + selfcal_library['Stop_Reason']=reason #selfcal_library[vis][solint]['Pass']=False selfcal_library[vis][solint]['Fail_Reason']=reason @@ -698,8 +697,8 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ 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-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] diff --git a/auto_selfcal/tests/test_auto_selfcal.py b/auto_selfcal/tests/test_auto_selfcal.py index b419c288..57b506c7 100644 --- a/auto_selfcal/tests/test_auto_selfcal.py +++ b/auto_selfcal/tests/test_auto_selfcal.py @@ -140,7 +140,7 @@ def test_on_github(tmp_path, request, zip_file, link): 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"]) + "am_smoothfactor", "Stop_Reason"]) assert difference_count == 0 diff --git a/auto_selfcal/weblog_creation.py b/auto_selfcal/weblog_creation.py index 1fa0f01b..c7a631e0 100644 --- a/auto_selfcal/weblog_creation.py +++ b/auto_selfcal/weblog_creation.py @@ -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][vis]])) or \ + (np.all(['Selfcal_Not_Attempted' in sclib[target][band][vis]['Stop_Reason'] for vis in sclib[target][band]['vislist']]))) and sclib[target][band]['final_solint']=='None': render_summary_table(htmlOut,sclib,target,band,directory=directory) continue htmlOut.writelines('Final Successful solint: '+str(sclib[target][band]['final_solint'])+'

\n') @@ -115,7 +117,7 @@ def generate_weblog(sclib,selfcal_plan,directory='weblog'): htmlOut.writelines('Noise Characteristics
\n') # Solint summary table - if 'Empty model' not in sclib[target][band]['Stop_Reason']: + if 'Empty model' not in sclib[target][band][sclib[target][band]['vislist'][0]]['Stop_Reason']: render_selfcal_solint_summary_table(htmlOut,sclib,target,band,selfcal_plan) # PER SPW STATS TABLE @@ -128,7 +130,7 @@ def generate_weblog(sclib,selfcal_plan,directory='weblog'): htmlOut.close() # Pages for each solint - if 'Empty model' not in sclib[target][band]['Stop_Reason']: + if 'Empty model' not in sclib[target][band][sclib[target][band]['vislist'][0]]['Stop_Reason']: render_per_solint_QA_pages(sclib,selfcal_plan,bands,directory=directory) @@ -566,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('\n') From 2643132bf19305240bc12208cf45f071d2793c51 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Wed, 13 May 2026 13:02:56 +0000 Subject: [PATCH 18/27] vis should be vislist --- auto_selfcal/weblog_creation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/auto_selfcal/weblog_creation.py b/auto_selfcal/weblog_creation.py index c7a631e0..61102b7b 100644 --- a/auto_selfcal/weblog_creation.py +++ b/auto_selfcal/weblog_creation.py @@ -88,7 +88,7 @@ def generate_weblog(sclib,selfcal_plan,directory='weblog'): 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][vis]])) or \ + 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 From 9eb76be714f158f41c3bacd17b374910255b6a43 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Tue, 2 Jun 2026 21:29:28 +0000 Subject: [PATCH 19/27] Remove print statements that are superseded by updates on this branch. --- auto_selfcal/run_selfcal.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/auto_selfcal/run_selfcal.py b/auto_selfcal/run_selfcal.py index 8a8f9812..fa2bf8f9 100644 --- a/auto_selfcal/run_selfcal.py +++ b/auto_selfcal/run_selfcal.py @@ -82,8 +82,6 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ for vis in selfcal_library['vislist']: print('selfcal_library[vis]["sub-fields-to-selfcal"] = ', selfcal_library[vis]['sub-fields-to-selfcal']) while iteration < len(selfcal_plan['solints']): - print(selfcal_plan['solints']) - print(selfcal_plan['solint_interval'][iteration]) 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:') From 39d1913da242a91e3b350d34c38fc99c63f67950 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Mon, 8 Jun 2026 20:58:13 +0000 Subject: [PATCH 20/27] Fix a few issues with the merge. --- auto_selfcal/prepare_cocal.py | 4 +--- auto_selfcal/run_selfcal.py | 11 ++++++----- auto_selfcal/tests/test_auto_selfcal.py | 1 - 3 files changed, 7 insertions(+), 9 deletions(-) diff --git a/auto_selfcal/prepare_cocal.py b/auto_selfcal/prepare_cocal.py index cf845663..88717500 100644 --- a/auto_selfcal/prepare_cocal.py +++ b/auto_selfcal/prepare_cocal.py @@ -55,11 +55,10 @@ def prepare_cocal(selfcal_library, selfcal_plan, inf_EB_gaincal_combine, inf_EB_ # Update the relevant lists if we are going to do a fallback mode. selfcal_plan[target][band]['solints'] += ["inf_EB_fb","inf_fb1","inf_fb2","inf_fb3"] selfcal_plan[target][band]['solmode'] += ["p","p","p","p"] - selfcal_plan[target][band]['solint_interval'] += ["inf","inf","inf","inf"] for vis in selfcal_library[target][band]['vislist']: - 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] = {} 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 @@ -89,7 +88,6 @@ def prepare_cocal(selfcal_library, selfcal_plan, inf_EB_gaincal_combine, inf_EB_ for cal_target in inf_fields[band]: selfcal_plan[target][band][vis]['solint_settings']["inf_fb3"]["preapply_gaintable_dict"][cal_target] = selfcal_plan[cal_target][band][vis.replace(sanitize_string(target), sanitize_string(cal_target))]['solint_settings']['inf_EB']['accepted_gaintable'] - selfcal_plan[target][band]['gaincal_combine'] += [selfcal_plan[target][band]['gaincal_combine'][0], selfcal_plan[target][band]['gaincal_combine'][1], selfcal_plan[target][band]['gaincal_combine'][1], selfcal_plan[target][band]['gaincal_combine'][1]] selfcal_plan[target][band]['applycal_mode'] += [selfcal_plan[target][band]['applycal_mode'][0], selfcal_plan[target][band]['applycal_mode'][1], selfcal_plan[target][band]['applycal_mode'][1], selfcal_plan[target][band]['applycal_mode'][1]] calibrators[band] = [inf_EB_fields[band], inf_fields[band], inf_fields[band], inf_fields[band]] diff --git a/auto_selfcal/run_selfcal.py b/auto_selfcal/run_selfcal.py index fa2bf8f9..e1c501c0 100644 --- a/auto_selfcal/run_selfcal.py +++ b/auto_selfcal/run_selfcal.py @@ -46,7 +46,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ if mode == "cocal": # Check whether there are suitable calibrators, otherwise skip this target/band. - include_targets, include_scans = triage_calibrators(vislist[0], target, band, calibrators[band][0]) + include_targets, include_scans = triage_calibrators(selfcal_library['vislist'][0], target, band, calibrators[band][0]) print('Co-calibrators: ',include_targets) print('Co-calibrator scans: ',include_scans) if include_targets == "": @@ -54,7 +54,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ selfcal_library['Stop_Reason'] += '; No suitable co-calibrators' return else: - for vis in vislist: + for vis in selfcal_library['vislist']: clearcal(vis=vis,addmodel=True) os.system('mv '+vis+' '+vis.replace('.ms','_orig.ms')) os.system('mv '+vis+'.flagversions '+vis.replace('.ms','_orig.ms.flagversions')) @@ -107,7 +107,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ calculate_inf_EB_fb_anyways = False preapply_targets_own_inf_EB = True # If there was no inf solint (e.g. because each source was observed only a single time, skip this as there are no gain tables to stick together. - if "inf" not in selfcal_plan['solints']: + if "inf" not in [selfcal_plan[vislist[0]]['solint_settings'][solint]['sub-name'] for solint in selfcal_plan['solints']]: iteration += 1 continue @@ -417,7 +417,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) or \ (((post_mosaic_RMS[fid] - mosaic_RMS[fid])/mosaic_RMS[fid] > 1.05 or \ (post_mosaic_RMS_NF[fid] - mosaic_RMS_NF[fid])/mosaic_RMS_NF[fid] > 1.05) and \ - selfcal_plan[fid]['solint_snr_per_field'][solint] > 5) + np.all([selfcal_plan[fid][vis]['solint_snr_per_field'][solint] > 5 for vis in selfcal_library['vislist-to-gaincal']])) if 'inf_EB' in solint or 'd' in solint: # If any of the fields succeed in the "strict" sense, then allow for minor reductions in the evaluation quantity in other @@ -498,7 +498,8 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ marginal_inf_EB_will_attempt_next_solint=False 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']]): + if np.all([selfcal_plan[vis]['solint_snr'][selfcal_plan['solints'][iteration+1]] < minsnr_to_proceed for vis in selfcal_library['vislist-to-gaincal']]) and \ + np.all([np.all([selfcal_plan[fid][vis]['solint_snr_per_field'][selfcal_plan['solints'][iteration+1]] < minsnr_to_proceed for fid in selfcal_library['sub-fields']]) for vis in selfcal_library['vislist-to-gaincal']]): marginal_inf_EB_will_attempt_next_solint = False else: marginal_inf_EB_will_attempt_next_solint = True diff --git a/auto_selfcal/tests/test_auto_selfcal.py b/auto_selfcal/tests/test_auto_selfcal.py index c9c06c7d..9a3d909b 100644 --- a/auto_selfcal/tests/test_auto_selfcal.py +++ b/auto_selfcal/tests/test_auto_selfcal.py @@ -228,7 +228,6 @@ def compare_two_dictionaries(dictionary1, dictionary2, path=[], exclude=[], tole else: alt_key = key - if type(dictionary1[alt_key]) == dict: difference_count += compare_two_dictionaries(dictionary1[alt_key], dictionary2[key], path.copy()+[key], exclude=exclude, tolerance=tolerance, key_map=key_map) else: From c1c5039342b888e0ba0020e2db53a06105f782b0 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Fri, 10 Jul 2026 07:24:56 -0400 Subject: [PATCH 21/27] Additional fixes to get the benchmark running through. --- auto_selfcal/prepare_selfcal.py | 2 +- auto_selfcal/run_selfcal.py | 7 +++++-- auto_selfcal/selfcal_helpers.py | 13 ++++++++++--- auto_selfcal/weblog_creation.py | 2 +- 4 files changed, 17 insertions(+), 7 deletions(-) diff --git a/auto_selfcal/prepare_selfcal.py b/auto_selfcal/prepare_selfcal.py index c53b66d2..54e3e5f2 100644 --- a/auto_selfcal/prepare_selfcal.py +++ b/auto_selfcal/prepare_selfcal.py @@ -164,7 +164,7 @@ 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) + selfcal_library[target][band][vis]['sub-fields-to-selfcal'].remove(fid) continue selfcal_library[target][band][fid][vis] = {} diff --git a/auto_selfcal/run_selfcal.py b/auto_selfcal/run_selfcal.py index e1c501c0..4ddedd0c 100644 --- a/auto_selfcal/run_selfcal.py +++ b/auto_selfcal/run_selfcal.py @@ -166,7 +166,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ ## if selfcal_library['final_solint'] != 'None': prev_solint = selfcal_library['final_solint'] - prev_iteration = selfcal_library[vislist[0]][prev_solint]['iteration'] + prev_iteration = selfcal_library[[vis for vis in selfcal_library['vislist-to-gaincal'] if prev_solint in selfcal_library[vis]][0]][prev_solint]['iteration'] nterms_changed = (len(glob.glob(sani_target+'_'+band+'_'+prev_solint+'_'+str(prev_iteration)+"_post.model.tt*")) < selfcal_library['nterms']) @@ -803,7 +803,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ print(vis+' now: ',selfcal_plan[vis]['solint_snr'][selfcal_plan['solints'][iteration+1]]) for fid in selfcal_library['sub-fields-to-selfcal']: - for vis in vislist: + for vis in selfcal_library[fid]['vislist']: if selfcal_plan['solints'][iteration+1] not in selfcal_plan[fid][vis]['solint_snr_per_field']: continue print('Field '+str(fid)+' '+vis+' was: ',selfcal_plan[fid][vis]['solint_snr_per_field'][selfcal_plan['solints'][iteration+1]]) @@ -857,6 +857,9 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ for vis in vislist: new_fields_to_selfcal = [] for fid in selfcal_library['sub-fields']: + if vis not in selfcal_library[fid]['vislist']: + continue + if mode == "cocal": if ("inf_EB" in selfcal_library[fid]['vislist'][0] and selfcal_library[fid][vis][coarsest_solint]["Pass"]) or selfcal_library[fid][selfcal_library[fid]['vislist'][0]][coarsest_solint+"_fb"]["Pass"]: new_fields_to_selfcal.append(fid) diff --git a/auto_selfcal/selfcal_helpers.py b/auto_selfcal/selfcal_helpers.py index 62312ac4..a7fb18e4 100644 --- a/auto_selfcal/selfcal_helpers.py +++ b/auto_selfcal/selfcal_helpers.py @@ -1817,6 +1817,7 @@ def get_SNR_self(selfcal_library,selfcal_plan,n_ant,inf_EB_gaincal_combine,inf_E for fid in selfcal_library[target][band]['sub-fields']: selfcal_plan[target][band][fid] = {} for vis in selfcal_library[target][band][fid]['vislist']: + print(target, band, vis, fid, selfcal_library[target][band][fid]['vislist']) selfcal_plan[target][band][fid][vis] = {} solints_per_vis = [solint for solint in selfcal_plan[target][band]['solints'] if solint in selfcal_plan[target][band][vis]['solint_settings']] selfcal_plan[target][band][fid][vis]['solint_snr_per_field'], selfcal_plan[target][band][fid][vis]['solint_snr_per_field_per_spw'], \ @@ -1893,15 +1894,21 @@ def get_SNR_self_individual(vislist,selfcal_library,n_ant,solints,solint_setting 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 + total_vis += 1 + if total_vis > 0: + mean_SNR_spw=mean_SNR_spw/total_vis + else: + mean_SNR_spw = 0.0 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(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) + if len(vislist) > 0: + mean_SNR_bb=mean_SNR_bb/len(vislist) + else: + mean_SNR_bb = 0.0 print('mean_SNR_bb',mean_SNR_bb,baseband) solint_snr_per_bb[solint][baseband]=mean_SNR_bb solint_snr[solint]=np.mean(SNR_self_EB) diff --git a/auto_selfcal/weblog_creation.py b/auto_selfcal/weblog_creation.py index 61102b7b..290d58f3 100644 --- a/auto_selfcal/weblog_creation.py +++ b/auto_selfcal/weblog_creation.py @@ -358,7 +358,7 @@ def render_selfcal_solint_summary_table(htmlOut,sclib,target,band,selfcal_plan): for quantity in ['Nsols_with_preflagged_data','Flagged_Sols_with_preflagged_data','Frac_Flagged_with_preflagged_data','Nsols_without_preflagged_data','Flagged_Sols_without_preflagged_data','Frac_Flagged_without_preflagged_data','Frac_sols_not_equiv_zero_per_spw','Frac_sols_not_equiv_zero_per_bb','SPW_Combine_Mode']: line='\n '+quantity+'\n' for solint in solint_list: - if solint in vis_keys and sclib[target][band][vis][solint]['Pass'] != 'None' and 'gaintable' in sclib[target][band][vis][solint]: + if solint in sclib[target][band][vis] and sclib[target][band][vis][solint]['Pass'] != 'None' and 'gaintable' in sclib[target][band][vis][solint]: # only evaluate last gaintable not the pre-apply table #gaintable=sclib[target][band][vis][solint]['gaintable'][len(sclib[target][band][vis][solint]['gaintable'])-1] #nflagged_sols, nsols=get_sols_flagged_solns(gaintable) From 3271b18e0c04ec905dc41cf2075c77db577fb76e Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Fri, 10 Jul 2026 07:28:12 -0400 Subject: [PATCH 22/27] Add solint_map for benchmark, along with excluding some relevant keys for getting the proper comparison across all of the changes on per-EB solints. Clean up printed output of compare_two_dictionaries, including adding a verbose option. --- auto_selfcal/tests/test_auto_selfcal.py | 50 ++++++++++++++++++------- 1 file changed, 36 insertions(+), 14 deletions(-) diff --git a/auto_selfcal/tests/test_auto_selfcal.py b/auto_selfcal/tests/test_auto_selfcal.py index 9a3d909b..26a1cc4e 100644 --- a/auto_selfcal/tests/test_auto_selfcal.py +++ b/auto_selfcal/tests/test_auto_selfcal.py @@ -36,7 +36,6 @@ pytest.param('disjoint_missing_spwfid_combo', id='disjoint_missing_spwfid_combo'), # tests spw-mapping, particularly in cases where spws and sub-fields are missing pytest.param('EE11.1.00032.S', id='EE11.1.00032.S'), # tests a dataset with spw names crossing a naming scheme change pytest.param('VLBA_CASA_Guide', id='VLBA_CASA_Guide'), # VLBA dataset - ] ) def test_benchmark(tmp_path, dataset): @@ -72,10 +71,28 @@ def test_benchmark(tmp_path, dataset): with open('selfcal_library.pickle', 'rb') as handle: selfcal_library2 = pickle.load(handle) - difference_count = compare_two_dictionaries(selfcal_library1, selfcal_library2, tolerance=1e-3,\ + with open('selfcal_plan.pickle', 'rb') as handle: + selfcal_plan = pickle.load(handle) + + solint_map = {} + for target in selfcal_library2: + for band in selfcal_library2[target]: + for vis in selfcal_library2[target][band]['vislist']: + for solint in selfcal_plan[target][band][vis]['solint_settings']: + if solint not in solint_map: + solint_map[solint] = [] + + mapped_solint = selfcal_plan[target][band][vis]['solint_settings'][solint]['sub-name'] + + solint_map[solint].append(mapped_solint) + print(solint_map) + + difference_count = compare_two_dictionaries(selfcal_library1, selfcal_library2, tolerance=1e-3, key_map=solint_map,\ exclude=['vislist_orig','field_str','imsize','flux_threshold','overlap_tol','bands_for_targets',\ 'am_dogrowprune','am_growiterations','am_lownoisethreshold','am_minbeamfrac',\ - 'am_noisethreshold','am_sidelobethreshold','am_smoothfactor','telescope']) + 'am_noisethreshold','am_sidelobethreshold','am_smoothfactor','telescope', "final_phase_solint", "final_solint", \ + "gaintable_final", "per_EB_SNR", "vislist-to-gaincal", "telescope", "gaintable","sub-fields-to-gaincal", \ + "sub-fields-to-selfcal"]) for msfile in starting_MS_files: print(f"rm -rf {os.path.basename(msfile)}") @@ -145,7 +162,7 @@ def test_on_github(tmp_path, request, zip_file, link): mapped_solint = selfcal_plan[target][band][vis]['solint_settings'][solint]['sub-name'] solint_map[solint].append(mapped_solint) - print(solint_map) + #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", @@ -175,7 +192,7 @@ def compare_values(list1, list2, tol=1e-3): else: return abs(list1 - list2) < abs(list1*tol) -def compare_two_dictionaries(dictionary1, dictionary2, path=[], exclude=[], tolerance=1e-3, key_map={}): +def compare_two_dictionaries(dictionary1, dictionary2, path=[], exclude=[], tolerance=1e-3, key_map={}, verbose=True): if isinstance(dictionary1, str): with open(dictionary1, 'rb') as handle: dictionary1 = pickle.load(handle) @@ -193,15 +210,17 @@ def compare_two_dictionaries(dictionary1, dictionary2, path=[], exclude=[], tole if key not in intersect_keys and key not in key_map and not np.any([key in key_map[k] for k in key_map]): if key not in dictionary1: - print('/'.join([str(p) for p in path])+"/"+key+" not in dictionary1") + if verbose: + print('/'.join([str(p) for p in path])+"/"+key+" not in dictionary1") else: - print('/'.join([str(p) for p in path])+"/"+key+" not in dictionary2") + if verbose: + print('/'.join([str(p) for p in path])+"/"+key+" not in dictionary2") difference_count += 1 continue elif key not in intersect_keys and key not in key_map and np.any([key in key_map[k] for k in key_map]): - print(f'key {key} has changed in dictionary2 and will be matched elsewhere') + #print(f'key {key} has changed in dictionary2 and will be matched elsewhere') continue try: @@ -211,25 +230,27 @@ def compare_two_dictionaries(dictionary1, dictionary2, path=[], exclude=[], tole continue if key in dictionary2 and not key in dictionary1 and key in key_map: - print(f'Checking whether key {key} has its name changed') + #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') + #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}") + pass + #print(f"Using alternative key {alt_key} to match with key {key}") else: - print(f"No match found in dictionary1, this is a difference") + if verbose: + print(f"No match for key {key} found in dictionary1, this is a difference") difference_count += 1 continue else: alt_key = key if type(dictionary1[alt_key]) == dict: - difference_count += compare_two_dictionaries(dictionary1[alt_key], dictionary2[key], path.copy()+[key], exclude=exclude, tolerance=tolerance, key_map=key_map) + difference_count += compare_two_dictionaries(dictionary1[alt_key], dictionary2[key], path.copy()+[key], exclude=exclude, tolerance=tolerance, key_map=key_map, verbose=verbose) else: value1 = np.array(dictionary1[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] @@ -243,7 +264,8 @@ def compare_two_dictionaries(dictionary1, dictionary2, path=[], exclude=[], tole value2.sort() if not compare_values(value1, value2, tol=tolerance): - print('/'.join([str(p) for p in path])+"/"+key, dictionary1[key], dictionary2[key]) + if verbose: + print('/'.join([str(p) for p in path])+"/"+key, dictionary1[key], dictionary2[key]) difference_count += 1 return difference_count From 64eef5c900d65ef2f11c028eeff79e00e827bbf1 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Wed, 22 Jul 2026 10:57:23 -0400 Subject: [PATCH 23/27] Properly handle marginal inf_EB cases for mosaics with the new per-EB solint framework. --- auto_selfcal/run_selfcal.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/auto_selfcal/run_selfcal.py b/auto_selfcal/run_selfcal.py index 4ddedd0c..6861a618 100644 --- a/auto_selfcal/run_selfcal.py +++ b/auto_selfcal/run_selfcal.py @@ -764,13 +764,17 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ # because the solint failed, because for mosaics we can keep trying the field as we clean deeper. If we set to False now, that wont happen. for fid in np.intersect1d(selfcal_library['sub-fields'],list(selfcal_library['sub-fields-fid_map'][vis].keys())): if (selfcal_library['final_solint'] == 'inf_EB' and selfcal_library['inf_EB_SNR_decrease']) or \ - ((not selfcal_library[vislist[0]][solint]['Pass'] or solint == 'int') and \ + (np.any([not selfcal_library[vis][solint]['Pass'] or selfcal_plan[vis]['solint_settings'][solint]['sub-name'] == 'int' for vis in selfcal_library[fid]['vislist-to-gaincal']]) and \ (selfcal_library[fid]['final_solint'] == 'inf_EB' and selfcal_library[fid]['inf_EB_SNR_decrease'])): - selfcal_library[fid]['SC_success']=False - selfcal_library[fid]['final_solint']='None' - for vis in selfcal_library[fid]['vislist']: - selfcal_library[fid][vis]['inf_EB']['Pass']=False # remove the success from inf_EB - selfcal_library[fid][vis]['inf_EB']['Fail_Reason']+=' with no successful solints later' # remove the success from inf_EB + + for vis in selfcal_library[fid]['vislist-to-gaincal']: + if not selfcal_library[vis][solint]['Pass'] or selfcal_plan[vis]['solint_settings'][solint]['sub-name'] == 'int': + selfcal_library[fid][vis]['inf_EB']['Pass']=False # remove the success from inf_EB + selfcal_library[fid][vis]['inf_EB']['Fail_Reason']+=' with no successful solints later' # remove the success from inf_EB + + if np.all([not selfcal_library[fid][vis]['inf_EB']['Pass'] for vis in selfcal_library[fid]['vislist']]): + selfcal_library[fid]['SC_success']=False + selfcal_library[fid]['final_solint']='None' for vis in vislist: applycal_wrapper(vis, target, band, solint, selfcal_library, From 6ffced76f655f65ce54efa9d570084a1396f9330 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Fri, 7 Aug 2026 08:18:15 -0400 Subject: [PATCH 24/27] Store Estimated_SNR_too_low reason on removal from vislist so it is handled properly in the weblog; Ensure sub-fields only have the gaintable applied if they are in sub-fields-to-selfcal for an EB. --- auto_selfcal/run_selfcal.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/auto_selfcal/run_selfcal.py b/auto_selfcal/run_selfcal.py index 6861a618..97a6b84d 100644 --- a/auto_selfcal/run_selfcal.py +++ b/auto_selfcal/run_selfcal.py @@ -130,6 +130,10 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ ' **************') remove_vis.append(vis) + selfcal_library[vis]['Stop_Reason']='Estimated_SNR_too_low_for_solint '+selfcal_plan['solints'][iteration] + for fid in selfcal_library['sub-fields-to-selfcal']: + selfcal_library[fid][vis]['Stop_Reason']='Estimated_SNR_too_low_for_solint '+selfcal_plan['solints'][iteration] + for rvis in remove_vis: vislist.remove(rvis) @@ -140,10 +144,6 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ print('****************Attempting amplitude selfcal*************') continue - 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 @@ -435,10 +435,10 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ print('****************Not all fields were successful, so re-applying and re-making _post image*************') for vis in vislist: applycal_wrapper(vis, target, band, solint, selfcal_library, - current=lambda f: f in field_by_field_success_dict and field_by_field_success_dict[f], - final=lambda f: (f not in field_by_field_success_dict or not field_by_field_success_dict[f]) and + current=lambda f: f in field_by_field_success_dict and field_by_field_success_dict[f] and f in selfcal_library[vis]['sub-fields-to-selfcal'], + final=lambda f: (f not in field_by_field_success_dict or not field_by_field_success_dict[f] or f not in selfcal_library[vis]['sub-fields-to-selfcal']) and selfcal_library[f]['SC_success'], - clear=lambda f: (f not in field_by_field_success_dict or not field_by_field_success_dict[f]) and + clear=lambda f: (f not in field_by_field_success_dict or not field_by_field_success_dict[f] or f not in selfcal_library[vis]['sub-fields-to-selfcal']) and not selfcal_library[f]['SC_success'], restore_flags='selfcal_starting_flags_'+sani_target+'_'+band) @@ -568,6 +568,8 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ selfcal_library[fid]['inf_EB_SNR_decrease']=False for vis in selfcal_library[fid]['vislist-to-gaincal']: + if fid not in selfcal_library[vis]['sub-fields-to-selfcal']: + continue # I think the below might not be correct for mosaics - it would set the gaintable even if fid is not in [vis]['sub-fields-to-selfcal']. # I think this needs: # if fid not in selfcal_library[vis]['sub-fields-to-selfcal']: @@ -700,7 +702,7 @@ def run_selfcal(selfcal_library, selfcal_plan, target, band, n_ants, \ reason=reason+'; ' reason=reason+'All sub-fields failed' for vis in vislist: - selfcal_library['Stop_Reason']=reason + selfcal_library[vis]['Stop_Reason']=reason #selfcal_library[vis][solint]['Pass']=False selfcal_library[vis][solint]['Fail_Reason']=reason From 2b1b2274892df8d605ff5bd89dfd4dab0da3736a Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Fri, 7 Aug 2026 08:18:54 -0400 Subject: [PATCH 25/27] Only store a solint_snr_per_spw if the spw is actually in an EB. --- auto_selfcal/selfcal_helpers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/auto_selfcal/selfcal_helpers.py b/auto_selfcal/selfcal_helpers.py index a7fb18e4..1c40721c 100644 --- a/auto_selfcal/selfcal_helpers.py +++ b/auto_selfcal/selfcal_helpers.py @@ -1897,9 +1897,9 @@ def get_SNR_self_individual(vislist,selfcal_library,n_ant,solints,solint_setting total_vis += 1 if total_vis > 0: mean_SNR_spw=mean_SNR_spw/total_vis + solint_snr_per_spw[solint][str(spw)]=mean_SNR_spw else: mean_SNR_spw = 0.0 - 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(vislist)): From 6e7feabf69243a40d2decf92d36c2c569acce813 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Fri, 7 Aug 2026 08:19:49 -0400 Subject: [PATCH 26/27] Show the result per EB, as these can differ, and as long as one EB has a Stop_Reason that is not None, show the stats. --- auto_selfcal/weblog_creation.py | 25 +++++++++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/auto_selfcal/weblog_creation.py b/auto_selfcal/weblog_creation.py index 290d58f3..7debcbb3 100644 --- a/auto_selfcal/weblog_creation.py +++ b/auto_selfcal/weblog_creation.py @@ -250,8 +250,26 @@ def render_selfcal_solint_summary_table(htmlOut,sclib,target,band,selfcal_plan): line += ' - \n' line += '\n' htmlOut.writelines(line) + htmlOut.writelines('\n Result: \n') + for vis in vislist: + line=f'\n {vis}: \n' + for solint in solint_list: + if solint in sclib[target][band][vis]: + if 'Pass' in sclib[target][band][vis][solint]: + if sclib[target][band][vis][solint]['Pass'] == False: + line+=' {} {}\n'.format('Fail',sclib[target][band][vis][solint]['Fail_Reason']) + elif sclib[target][band][vis][solint]['Pass'] == 'None': + line+=' {} {}\n'.format('Not attempted',sclib[target][band][vis][solint]['Fail_Reason']) + else: + line+=' {}\n'.format('Pass') + else: + line+=' {}\n'.format('None') + else: + line += ' - \n' + line += '\n' + htmlOut.writelines(line) htmlOut.writelines('\n Selfcal stats: \n') - quantities=['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'] + quantities=['intflux_final','intflux_improvement','SNR_final','SNR_Improvement','SNR_NF_final','SNR_NF_Improvement','RMS_final','RMS_Improvement','RMS_NF_final','RMS_NF_Improvement','Beam_Ratio','clean_threshold','Plots'] for key in quantities: if key =='Pass': line='\n Result: \n' @@ -283,7 +301,10 @@ def render_selfcal_solint_summary_table(htmlOut,sclib,target,band,selfcal_plan): 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] + if np.any([sclib[target][band][vis][solint]['Pass'] != 'None' for vis in vislist]): + ivis = np.where([solint in sclib[target][band][vis] and sclib[target][band][vis][solint]['Pass'] != 'None' for vis in vislist])[0][0] + else: + ivis = np.where([solint in sclib[target][band][vis] for vis in vislist])[0][0] else: ivis = len(vislist)-1 From 2a37316fea02712ffbedb2dfeff901e9c6f92926 Mon Sep 17 00:00:00 2001 From: Patrick Sheehan Date: Fri, 7 Aug 2026 08:22:57 -0400 Subject: [PATCH 27/27] Use uniform_solints=True for tests at the moment, and exclude a few additional selfcal_library entries that are/can be handled differently even if the outcome is the same. --- auto_selfcal/tests/test_auto_selfcal.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/auto_selfcal/tests/test_auto_selfcal.py b/auto_selfcal/tests/test_auto_selfcal.py index 26a1cc4e..f4737dc1 100644 --- a/auto_selfcal/tests/test_auto_selfcal.py +++ b/auto_selfcal/tests/test_auto_selfcal.py @@ -61,7 +61,7 @@ def test_benchmark(tmp_path, dataset): if 'VLBA' in dataset: job = ex.submit(auto_selfcal, sort_targets_and_EBs=True, weblog=True, parallel=True, iscalibrator=True, targets='J1154+6022',applytargets='J1203+6031',imsize=640,cell='0.0002arcsec') else: - job = ex.submit(auto_selfcal, sort_targets_and_EBs=True, weblog=True, parallel=True) + job = ex.submit(auto_selfcal, sort_targets_and_EBs=True, weblog=True, parallel=True, uniform_solints=True) job.wait() assert job.state in ['DONE','COMPLETED'] @@ -91,8 +91,8 @@ def test_benchmark(tmp_path, dataset): 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', "final_phase_solint", "final_solint", \ - "gaintable_final", "per_EB_SNR", "vislist-to-gaincal", "telescope", "gaintable","sub-fields-to-gaincal", \ - "sub-fields-to-selfcal"]) + "gaintable_final", "per_EB_SNR", "per_scan_SNR", "vislist-to-gaincal", "telescope", "gaintable","sub-fields-to-gaincal", \ + "sub-fields-to-selfcal", "Stop_Reason"]) for msfile in starting_MS_files: print(f"rm -rf {os.path.basename(msfile)}")