From 7a2f98a6e7725b91fb7b289c973549da8793488b Mon Sep 17 00:00:00 2001 From: Nichole Barry Date: Wed, 8 Jan 2025 15:22:51 +1100 Subject: [PATCH 1/7] bandpass to metadata --- fhd_core/calibration/fhd_struct_init_cal.pro | 2 +- fhd_core/calibration/vis_calibrate.pro | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/fhd_core/calibration/fhd_struct_init_cal.pro b/fhd_core/calibration/fhd_struct_init_cal.pro index 4953a967..7fb0e455 100644 --- a/fhd_core/calibration/fhd_struct_init_cal.pro +++ b/fhd_core/calibration/fhd_struct_init_cal.pro @@ -71,7 +71,7 @@ auto_scale=Fltarr(2) cal_struct={n_pol:n_pol, n_freq:n_freq, n_tile:n_tile, n_time:n_time, uu:u_loc, vv:v_loc, $ auto_initialize:auto_initialize, max_iter:max_cal_iter, phase_iter:phase_fit_iter, $ tile_A:tile_A, tile_B:tile_B, tile_names:tile_names, bin_offset:bin_offset, freq:freq, gain:gain_arr_ptr, $ - adaptive_gain:use_adaptive_calibration_gain, base_gain:calibration_base_gain, $ + gain_bandpass:Pointer_copy(gain_arr_ptr), adaptive_gain:use_adaptive_calibration_gain, base_gain:calibration_base_gain, $ gain_residual:gain_residual, auto_scale:auto_scale, auto_params:auto_params, cross_phase:0.0, stokes_mix_phase:0.0, $ min_cal_baseline:min_cal_baseline, max_cal_baseline:max_cal_baseline, n_vis_cal:n_vis_cal, $ time_avg:cal_time_average, min_solns:min_cal_solutions, ref_antenna:ref_antenna, $ diff --git a/fhd_core/calibration/vis_calibrate.pro b/fhd_core/calibration/vis_calibrate.pro index 4b3259ed..f9694f60 100644 --- a/fhd_core/calibration/vis_calibrate.pro +++ b/fhd_core/calibration/vis_calibrate.pro @@ -206,6 +206,7 @@ FUNCTION vis_calibrate,vis_ptr,cal,obs,status_str,psf,params,jones,vis_weight_pt cal=vis_cal_combine(cal_polyfit,cal_bandpass) ENDIF ELSE cal=cal_bandpass + cal.gain_bandpass = cal_bandpass.gain IF Keyword_Set(auto_ratio_calibration) THEN BEGIN cal=cal_auto_ratio(obs,cal,auto_ratio=auto_ratio,auto_tile_i=auto_tile_i,/remultiply) From ee9f50fd29d1a45d2b63db52f31a74e647aff98c Mon Sep 17 00:00:00 2001 From: Nichole Barry Date: Wed, 8 Jan 2025 15:26:42 +1100 Subject: [PATCH 2/7] handle nans in cal plots --- fhd_output/calibration_plots/plot_cals.pro | 6 +++--- fhd_output/calibration_plots/plot_cals_sub.pro | 10 +++++----- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/fhd_output/calibration_plots/plot_cals.pro b/fhd_output/calibration_plots/plot_cals.pro index 967544c9..08123a8e 100644 --- a/fhd_output/calibration_plots/plot_cals.pro +++ b/fhd_output/calibration_plots/plot_cals.pro @@ -110,9 +110,9 @@ for file_i=0, page_num-1 do begin if page_num GT 1 then gains1_res=gains1_res[*,file_i*page_tile_max:(file_i+1)*page_tile_max<(size(gains1_res))[2]-1] gains1_orig=gains1_res+gains1 gains1_res_abs=Abs(gains1_orig)-Abs(gains1) - max_amp_res = mean(abs([gains0_res_abs,gains1_res_abs]),/NAN) + 3*stddev(abs([gains0_res_abs,gains1_res_abs]),/NAN) - ENDIF ELSE max_amp_res = mean(abs(gains0_res_abs),/NAN) + 3*stddev(abs(gains0_res_abs),/NAN) - + max_amp_res = mean(abs([gains0_res_abs,gains1_res_abs]),/nan) + 3*stddev(abs([gains0_res_abs,gains1_res_abs]),/nan) + ENDIF ELSE max_amp_res = mean(abs(gains0_res_abs),/nan) + 3*stddev(abs(gains0_res_abs),/nan) + ;plot amplitude residuals IF max_amp_res GT 0 THEN BEGIN plot_cals_sub,freq,gains0_orig,gains1_orig,filename=orig_amp_filename,tile_A=tile_A,tile_B=tile_B,tile_use=tile_use,tile_exist=tile_exist,tile_names=tile_names,$ diff --git a/fhd_output/calibration_plots/plot_cals_sub.pro b/fhd_output/calibration_plots/plot_cals_sub.pro index cd7ba840..1a794b4a 100644 --- a/fhd_output/calibration_plots/plot_cals_sub.pro +++ b/fhd_output/calibration_plots/plot_cals_sub.pro @@ -23,14 +23,14 @@ IF Keyword_Set(phase) THEN BEGIN ENDIF ELSE BEGIN IF N_Elements(yrange) NE 2 THEN BEGIN IF n_pol GT 1 THEN BEGIN - yrange=Minmax(abs([gains_A[*,tile_use_i],gains_B[*,tile_use_i]])) + yrange=Minmax(abs([gains_A[*,tile_use_i],gains_B[*,tile_use_i]]),/nan) ENDIF ELSE BEGIN - yrange=Minmax(abs(gains_A[*,tile_use_i])) + yrange=Minmax(abs(gains_A[*,tile_use_i]),/nan) ENDELSE IF (gain_type LE 5) OR Keyword_Set(real_vs_imaginary) THEN BEGIN - IF n_pol GT 1 THEN max_amp = mean(abs([gains_A[*,tile_use_i],gains_B[*,tile_use_i]]),/NAN) + 5*stddev(abs([gains_A[*,tile_use_i],gains_B[*,tile_use_i]]),/NAN) $ - ELSE max_amp = Mean(abs(gains_A[*,tile_use_i]),/NAN) + 5*stddev(abs(gains_A[*,tile_use_i]),/NAN) + IF n_pol GT 1 THEN max_amp = mean(abs([gains_A[*,tile_use_i],gains_B[*,tile_use_i]]),/nan) + 5*stddev(abs([gains_A[*,tile_use_i],gains_B[*,tile_use_i]]),/nan) $ + ELSE max_amp = Mean(abs(gains_A[*,tile_use_i]),/nan) + 5*stddev(abs(gains_A[*,tile_use_i]),/nan) amp_range=Floor(Alog10(max_amp)) IF amp_range LT 0 THEN amp_digits=1. ELSE amp_digits=2. amp_test=max_amp/10.^(amp_range-amp_digits) @@ -104,4 +104,4 @@ ENDFOR cgtext,.4,max(plot_pos[*,3]+height/4),obsname,/normal cgPS_Close,/png,Density=300,Resize=cal_plot_resize,/allow_transparent,/nomessage -END +END \ No newline at end of file From 31ac79fe547230c2802b8886c1b55ba5d4f57ed2 Mon Sep 17 00:00:00 2001 From: Nichole Barry Date: Thu, 9 Jan 2025 14:24:37 +1100 Subject: [PATCH 3/7] model transfer generalization --- fhd_core/calibration/vis_calibrate.pro | 2 +- .../source_modeling/vis_model_transfer.pro | 75 ++++++++++++++++++- fhd_core/source_modeling/vis_source_model.pro | 2 +- 3 files changed, 76 insertions(+), 3 deletions(-) diff --git a/fhd_core/calibration/vis_calibrate.pro b/fhd_core/calibration/vis_calibrate.pro index f9694f60..353b3da2 100644 --- a/fhd_core/calibration/vis_calibrate.pro +++ b/fhd_core/calibration/vis_calibrate.pro @@ -96,7 +96,7 @@ FUNCTION vis_calibrate,vis_ptr,cal,obs,status_str,psf,params,jones,vis_weight_pt ;; Option to transfer pre-made and unflagged model visbilities if keyword_set(model_transfer) then begin - vis_model_arr = vis_model_transfer(obs,model_transfer) + vis_model_arr = vis_model_transfer(obs,params,model_transfer) endif RETURN,vis_cal diff --git a/fhd_core/source_modeling/vis_model_transfer.pro b/fhd_core/source_modeling/vis_model_transfer.pro index 3f2aedf7..07710be1 100644 --- a/fhd_core/source_modeling/vis_model_transfer.pro +++ b/fhd_core/source_modeling/vis_model_transfer.pro @@ -1,4 +1,7 @@ -function vis_model_transfer,obs,model_transfer +function vis_model_transfer,obs,params,model_transfer + + data_nbaselines = obs.nbaselines + data_n_time = obs.n_time ;; Option to transfer pre-made and unflagged model visbilities vis_model_arr=PTRARR(obs.n_pol,/allocate) @@ -11,6 +14,76 @@ function vis_model_transfer,obs,model_transfer print, "Model visibilities transferred from " + transfer_name endfor + ;; Get the params file associated with the model visibilities + if file_test(model_transfer+'/'+obs.obsname+'_params.sav') then begin + params_model = getvar_savefile(model_transfer+'/'+obs.obsname+'_params.sav','params') + endif else begin + if file_test(file_dirname(model_transfer)+'/metadata/'+obs.obsname+'_params.sav') then begin + params_model = getvar_savefile(file_dirname(model_transfer)+'/metadata/'+obs.obsname+'_params.sav','params') + endif else message, 'No params file found in model transfer directory.' + endelse + model_n_time = n_elements(uniq(params_model.time)) + model_nbaselines = (uniq(params_model.time))[0]+1 + + if model_n_time NE data_n_time then begin + ;; Exclude flagged times from the model visibilities if the number of time steps do not match + + data_Jdate = (*obs.baseline_info).Jdate + time_res_Jdate = (double(obs.time_res)/2.) / (24.*3600.) ;offset to the center of a time step + model_time = params_model.time[uniq(params_model.time)] + data_time = params.time[uniq(params.time)] + tolerance = 1E-5 + + ; Match times betweeen model and data + ; Option 1: You are working with a model made by FHD. Thus the timing convention in the params is the same + ; Option 2: You are working with a model made by WODEN. Thus the timing convention in the params is different + ; and you need to convert the model Jdate to the data Jdate within precision. + + ; Option 1 + matched_times_opt1 = intarr(data_n_time) + for data_time_i=0, data_n_time-1 do begin + min_val= min(abs(model_time - data_time[data_time_i]), min_ind) + if min_val LT tolerance then matched_times_opt1[data_time_i] = min_ind + endfor + + ; Option 2 + matched_times_opt2 = intarr(data_n_time) + for data_time_i=0, data_n_time-1 do begin + min_val = min(abs(model_time - (data_Jdate[data_time_i]-time_res_Jdate)), min_ind) + if min_val LT tolerance then matched_times_opt2[data_time_i] = min_ind + endfor + + ; Find which option was more successful in matching the times + if total(matched_times_opt1 < 1) GT total(matched_times_opt2 < 1) then matched_times = matched_times_opt1 + if total(matched_times_opt1 < 1) LT total(matched_times_opt2 < 1) then matched_times = matched_times_opt2 + + endif else matched_times = indgen(model_n_time) + + ;An arbitrary larger number than tiles to make an index array + baseline_mod = ULONG(obs.n_tile*2) + ;Create a baseline array. Could use params.baseline_arr, but there is no guarentee of model gen origin + data_baseline_index = params.antenna1*baseline_mod + params.antenna2 + model_baseline_index = params_model.antenna1*baseline_mod + params_model.antenna2 + + ;Initialize matched model pointer array + matched_model = PTRARR(obs.n_pol,/allocate) + for pol_i=0, obs.n_pol-1 do begin + if size(*vis_model_arr[pol_i], /type) eq 6 then *matched_model[pol_i] = complex(FLTARR(obs.n_freq, data_nbaselines)) + if size(*vis_model_arr[pol_i], /type) eq 9 then *matched_model[pol_i] = dcomplex(DBLARR(obs.n_freq, data_nbaselines)) + endfor + + ; In each matched timestep, match the baselines and fill a new, matched model array + for time_i=0, data_n_time-1 do begin + ;suba is the subset of indices in the first array that are also in the second. + ;subb is the subset of indices in the second array that are also in the first. + match, data_baseline_index[time_i*data_nbaselines:(time_i+1)*data_nbaselines-1], $ + model_baseline_index[matched_times[time_i]*model_nbaselines:(matched_times[time_i]+1)*model_nbaselines-1], $ + suba, subb + for pol_i=0, obs.n_pol-1 do (*matched_model[pol_i])[*,suba+time_i*(data_nbaselines)] = (*vis_model_arr[pol_i])[*,subb+matched_times[time_i]*model_nbaselines] + endfor + + for pol_i=0, obs.n_pol-1 do vis_model_arr[pol_i] = Pointer_copy(matched_model[pol_i]) + return, vis_model_arr end diff --git a/fhd_core/source_modeling/vis_source_model.pro b/fhd_core/source_modeling/vis_source_model.pro index 4af865b2..7e51c1b8 100644 --- a/fhd_core/source_modeling/vis_source_model.pro +++ b/fhd_core/source_modeling/vis_source_model.pro @@ -6,7 +6,7 @@ FUNCTION vis_source_model,skymodel, obs, status_str, psf, params, vis_weight_ptr ;; Option to transfer pre-made and unflagged model visbilities if keyword_set(model_transfer) then begin - vis_arr = vis_model_transfer(obs,model_transfer) + vis_arr = vis_model_transfer(obs,params,model_transfer) return, vis_arr endif From 4c6911088c51bbf765fd8c2bfda37b83dc22bf49 Mon Sep 17 00:00:00 2001 From: Nichole Barry Date: Thu, 9 Jan 2025 17:05:12 +1100 Subject: [PATCH 4/7] array size fix --- fhd_core/source_modeling/vis_model_transfer.pro | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fhd_core/source_modeling/vis_model_transfer.pro b/fhd_core/source_modeling/vis_model_transfer.pro index 07710be1..4b28a2e7 100644 --- a/fhd_core/source_modeling/vis_model_transfer.pro +++ b/fhd_core/source_modeling/vis_model_transfer.pro @@ -68,8 +68,8 @@ function vis_model_transfer,obs,params,model_transfer ;Initialize matched model pointer array matched_model = PTRARR(obs.n_pol,/allocate) for pol_i=0, obs.n_pol-1 do begin - if size(*vis_model_arr[pol_i], /type) eq 6 then *matched_model[pol_i] = complex(FLTARR(obs.n_freq, data_nbaselines)) - if size(*vis_model_arr[pol_i], /type) eq 9 then *matched_model[pol_i] = dcomplex(DBLARR(obs.n_freq, data_nbaselines)) + if size(*vis_model_arr[pol_i], /type) eq 6 then *matched_model[pol_i] = complex(FLTARR(obs.n_freq, data_nbaselines * data_n_time)) + if size(*vis_model_arr[pol_i], /type) eq 9 then *matched_model[pol_i] = dcomplex(DBLARR(obs.n_freq, data_nbaselines * data_n_time)) endfor ; In each matched timestep, match the baselines and fill a new, matched model array From 10a26d24ffb40a165cb21f2ce6fd917173c32489 Mon Sep 17 00:00:00 2001 From: Nichole Barry Date: Wed, 22 Jan 2025 09:20:11 +1100 Subject: [PATCH 5/7] no_png option in bandpass plots --- fhd_core/calibration/vis_cal_bandpass.pro | 5 +++-- fhd_core/calibration/vis_calibrate.pro | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/fhd_core/calibration/vis_cal_bandpass.pro b/fhd_core/calibration/vis_cal_bandpass.pro index bd9f4703..e4d44d26 100644 --- a/fhd_core/calibration/vis_cal_bandpass.pro +++ b/fhd_core/calibration/vis_cal_bandpass.pro @@ -1,6 +1,7 @@ FUNCTION vis_cal_bandpass,cal,obs,params,cal_remainder=cal_remainder,file_path_fhd=file_path_fhd,cable_bandpass_fit=cable_bandpass_fit,$ bandpass_directory=bandpass_directory,tile_use=tile_use,calibration_bandpass_cable_exclude=calibration_bandpass_cable_exclude,$ - cal_bp_transfer=cal_bp_transfer,uvfits_version=uvfits_version,uvfits_subversion=uvfits_subversion,auto_ratio_calibration=auto_ratio_calibration,Extra=extra + cal_bp_transfer=cal_bp_transfer,uvfits_version=uvfits_version,uvfits_subversion=uvfits_subversion,auto_ratio_calibration=auto_ratio_calibration,$ + no_png=no_png,Extra=extra ;This function is version 1 of calibrating each group of tiles with similar cable lengths per observation. ;Extract needed elements from the input structures @@ -177,7 +178,7 @@ FUNCTION vis_cal_bandpass,cal,obs,params,cal_remainder=cal_remainder,file_path_f ENDELSE - IF Keyword_Set(file_path_fhd) THEN bandpass_plots,obs,bandpass_arr,file_path_fhd=file_path_fhd,$ + IF Keyword_Set(file_path_fhd) AND ~Keyword_Set(no_png) THEN bandpass_plots,obs,bandpass_arr,file_path_fhd=file_path_fhd,$ cable_length_ref=cable_length_ref,tile_use_arr=tile_use_arr RETURN,cal_bandpass diff --git a/fhd_core/calibration/vis_calibrate.pro b/fhd_core/calibration/vis_calibrate.pro index 353b3da2..93113564 100644 --- a/fhd_core/calibration/vis_calibrate.pro +++ b/fhd_core/calibration/vis_calibrate.pro @@ -199,7 +199,7 @@ FUNCTION vis_calibrate,vis_ptr,cal,obs,status_str,psf,params,jones,vis_weight_pt ENDIF cal_bandpass=vis_cal_bandpass(cal,obs,params,cal_remainder=cal_remainder,auto_ratio_calibration=auto_ratio_calibration,$ - file_path_fhd=file_path_fhd,_Extra=extra) + file_path_fhd=file_path_fhd,no_png=no_png,_Extra=extra) IF Keyword_Set(calibration_polyfit) THEN BEGIN cal_polyfit=vis_cal_polyfit(cal_remainder,obs,auto_ratio=auto_ratio,_Extra=extra) From caf8c551a457191dd665e55b0c6f96bd1a991ea9 Mon Sep 17 00:00:00 2001 From: Nichole Barry Date: Wed, 12 Feb 2025 10:12:57 +1100 Subject: [PATCH 6/7] add error messages --- fhd_core/source_modeling/vis_model_transfer.pro | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/fhd_core/source_modeling/vis_model_transfer.pro b/fhd_core/source_modeling/vis_model_transfer.pro index 4b28a2e7..9496dd72 100644 --- a/fhd_core/source_modeling/vis_model_transfer.pro +++ b/fhd_core/source_modeling/vis_model_transfer.pro @@ -54,10 +54,15 @@ function vis_model_transfer,obs,params,model_transfer endfor ; Find which option was more successful in matching the times + ; Error if no matching times found or if matching times does not obviously cover all data times + if (total(matched_times_opt1 < 1) EQ 0) and (total(matched_times_opt2 < 1) EQ 0) then $ + message, 'ERROR: No matching times found between transferred model and data.' if total(matched_times_opt1 < 1) GT total(matched_times_opt2 < 1) then matched_times = matched_times_opt1 if total(matched_times_opt1 < 1) LT total(matched_times_opt2 < 1) then matched_times = matched_times_opt2 + if total(matched_times < 1) NE data_n_time then $ + message, 'ERROR: Number of matching times found between transferred model and data does not equal the number of timesteps in data.' - endif else matched_times = indgen(model_n_time) + endif else matched_times = indgen(model_n_time) ;An arbitrary larger number than tiles to make an index array baseline_mod = ULONG(obs.n_tile*2) @@ -73,6 +78,7 @@ function vis_model_transfer,obs,params,model_transfer endfor ; In each matched timestep, match the baselines and fill a new, matched model array + ; Unmatched times in the transferred model are not included in the matched model array for time_i=0, data_n_time-1 do begin ;suba is the subset of indices in the first array that are also in the second. ;subb is the subset of indices in the second array that are also in the first. From 5eef80fb623a970fb09b443d8689ae2c246df8da Mon Sep 17 00:00:00 2001 From: Nichole Barry Date: Thu, 20 Feb 2025 16:12:07 +1100 Subject: [PATCH 7/7] add option 3 --- .../source_modeling/vis_model_transfer.pro | 40 ++++++++++++++----- 1 file changed, 30 insertions(+), 10 deletions(-) diff --git a/fhd_core/source_modeling/vis_model_transfer.pro b/fhd_core/source_modeling/vis_model_transfer.pro index 9496dd72..d590d2c0 100644 --- a/fhd_core/source_modeling/vis_model_transfer.pro +++ b/fhd_core/source_modeling/vis_model_transfer.pro @@ -36,30 +36,50 @@ function vis_model_transfer,obs,params,model_transfer ; Match times betweeen model and data ; Option 1: You are working with a model made by FHD. Thus the timing convention in the params is the same - ; Option 2: You are working with a model made by WODEN. Thus the timing convention in the params is different - ; and you need to convert the model Jdate to the data Jdate within precision. + ; Option 2: You are working with a model made by WODEN above version 1.4. This is different from FHD convention + ; by half a timestep. + ; Option 3: You are working with a model made by WODEN below version 1.4. Thus the timing convention in the + ; params is different and you need to convert the model Jdate to the data Jdate within precision. ; Option 1 - matched_times_opt1 = intarr(data_n_time) + matched_times_opt1 = intarr(data_n_time)-1 for data_time_i=0, data_n_time-1 do begin min_val= min(abs(model_time - data_time[data_time_i]), min_ind) if min_val LT tolerance then matched_times_opt1[data_time_i] = min_ind endfor - + ; Option 2 - matched_times_opt2 = intarr(data_n_time) + matched_times_opt2 = intarr(data_n_time)-1 for data_time_i=0, data_n_time-1 do begin - min_val = min(abs(model_time - (data_Jdate[data_time_i]-time_res_Jdate)), min_ind) + min_val= min(abs(model_time - (data_time[data_time_i]-time_res_Jdate)), min_ind) if min_val LT tolerance then matched_times_opt2[data_time_i] = min_ind endfor + ; Option 3 + matched_times_opt3 = intarr(data_n_time)-1 + for data_time_i=0, data_n_time-1 do begin + min_val = min(abs(model_time - (data_Jdate[data_time_i]-time_res_Jdate)), min_ind) + if min_val LT tolerance then matched_times_opt3[data_time_i] = min_ind + endfor + ; Find which option was more successful in matching the times + ; Initialization of matched_times arrays above to -1 allows up to check if all data times are matched (including the 0th index) ; Error if no matching times found or if matching times does not obviously cover all data times - if (total(matched_times_opt1 < 1) EQ 0) and (total(matched_times_opt2 < 1) EQ 0) then $ + temp = where(matched_times_opt1 NE -1, n_matched_opt1) + temp = where(matched_times_opt2 NE -1, n_matched_opt2) + temp = where(matched_times_opt3 NE -1, n_matched_opt3) + if (n_matched_opt1 EQ 0) and (n_matched_opt2 EQ 0) and (n_matched_opt3 EQ 0) then $ message, 'ERROR: No matching times found between transferred model and data.' - if total(matched_times_opt1 < 1) GT total(matched_times_opt2 < 1) then matched_times = matched_times_opt1 - if total(matched_times_opt1 < 1) LT total(matched_times_opt2 < 1) then matched_times = matched_times_opt2 - if total(matched_times < 1) NE data_n_time then $ + + max_matched = max([n_matched_opt1, n_matched_opt2, n_matched_opt3],max_matched_ind) + CASE max_matched_ind OF + 0: matched_times = matched_times_opt1 + 1: matched_times = matched_times_opt2 + 2: matched_times = matched_times_opt3 + ENDCASE + + temp = where(matched_times NE -1, n_matched) + if n_matched NE data_n_time then $ message, 'ERROR: Number of matching times found between transferred model and data does not equal the number of timesteps in data.' endif else matched_times = indgen(model_n_time)