Add redundant claibration option using visibility covariances

fhd_core/calibration/apply_redundant_cal_correction.pro

@@ -0,0 +1,20 @@
+FUNCTION apply_redundant_cal_correction, vis_model, freq_i=freq_i, redundant_cal_correction=redundant_cal_correction,$
+    baseline_inds=baseline_inds, use_conjugate=use_conjugate, redundant_calibration_weight=redundant_calibration_weight
+
+    IF N_Elements(redundant_calibration_weight) EQ 0 THEN redundant_calibration_weight=1.
+
+    IF N_Elements(freq_i) GT 0 THEN BEGIN
+        redundant_cal_correction_use = Reform(redundant_cal_correction[freq_i, *])
+        IF N_Elements(baseline_inds) GT 0 THEN $
+            redundant_cal_correction_use = redundant_cal_correction_use[baseline_inds]
+        IF Keyword_Set(use_conjugate) THEN redundant_cal_correction_use = [redundant_cal_correction_use, Conj(redundant_cal_correction_use)]
+    ENDIF ELSE BEGIN
+        IF N_Elements(baseline_inds) GT 0 THEN $
+            redundant_cal_correction_use = redundant_cal_correction[*, baseline_inds] ELSE $
+            redundant_cal_correction_use = redundant_cal_correction
+    ENDELSE
+
+    vis_model_corrected = vis_model + redundant_cal_correction_use*redundant_calibration_weight
+
+    RETURN, vis_model_corrected
+END

fhd_core/calibration/calculate_baseline_covariance.pro

@@ -0,0 +1,128 @@
+FUNCTION calculate_baseline_covariance,obs, psf, cal, freq_i, pol_i, baseline_inds=baseline_inds
+    ; Calculate the fraction of power shared between every baseline pair
+    n_pol = cal.n_pol
+    frequency = (*obs.baseline_info).freq[freq_i]
+
+    psf_dim=psf.dim
+    psf_resolution=psf.resolution
+    beam_arr=*psf.beam_ptr
+    freq_bin_i=(*obs.baseline_info).fbin_i
+    fbin_use=freq_bin_i[freq_i]
+    kbinsize=obs.kpix
+    kx_arr=cal.uu[baseline_inds]/kbinsize
+    ky_arr=cal.vv[baseline_inds]/kbinsize
+    conj_i=where(ky_arr GT 0,n_conj)
+    conj_flag=intarr(N_Elements(ky_arr)) 
+    IF n_conj GT 0 THEN BEGIN
+        conj_flag[conj_i]=1
+        kx_arr[conj_i]=-kx_arr[conj_i]
+        ky_arr[conj_i]=-ky_arr[conj_i]
+    ENDIF
+;    ; As with the rest of calibration, append the complex conjugate at the mirror location
+;    kx_arr = [kx_arr, -kx_arr]
+;    ky_arr = [ky_arr, -ky_arr]
+;    baseline_inds2 = [baseline_inds, baseline_inds]
+    xcen=frequency*Temporary(kx_arr)
+    ycen=frequency*Temporary(ky_arr)
+    n_baselines = N_Elements(baseline_inds)
+;    n_baselines2 = N_Elements(baseline_inds2)
+
+    psf_dim2 = 2*psf_dim + 1
+    init_arr=Complexarr(psf_dim2,psf_dim2)
+    init_arr[Ceil(psf_dim/2):Ceil(psf_dim/2)+psf_dim-1, Ceil(psf_dim/2):Ceil(psf_dim/2)+psf_dim-1] = 1
+    psf_base_inds = where(init_arr)
+    arr_type=Size(init_arr,/type)
+    psf_dim3 = psf_dim2*psf_dim2
+    ; Sparse matrix format, see holo_mapfn_convert.pro
+    sa=Ptrarr(n_baselines)
+    ija=Ptrarr(n_baselines)
+    overlap_test = Fltarr(psf_dim)
+    covariance_box = Make_array(psf_dim3, type=arr_type)
+    covariance_box[psf_base_inds] = *(*beam_arr[pol_i,fbin_use,0])[0,0]
+    covariance_box /= Total(Abs(covariance_box)^2.)
+    FOR pix=0,psf_dim-1 DO BEGIN
+        ind_offset = pix +pix*psf_dim2
+        covariance_box2 = Make_array(psf_dim3, type=arr_type)
+        covariance_box2[psf_base_inds + ind_offset] = Conj(*(*beam_arr[pol_i,fbin_use,0])[0,0])
+        overlap_test[pix] = Total(Abs(covariance_box*covariance_box2))
+    ENDFOR
+    overlap_i = where(overlap_test GE cal.covariance_threshold, n_overlap)
+    covariant_bin = n_overlap - 1
+    ; Pad the coordinates so that there is always a buffer of one bin around the edges
+    x_min = Min(xcen) - covariant_bin
+    y_min = Min(ycen) - covariant_bin
+    x_max = Max(xcen) + covariant_bin
+    y_max = Max(ycen) + covariant_bin
+    xsize = Ceil(x_max - x_min)
+    ysize = Ceil(y_max - y_min)
+    xsize = covariant_bin*Ceil(xsize/covariant_bin)
+    ysize = covariant_bin*Ceil(ysize/covariant_bin)
+    lookup_dim = xsize/covariant_bin
+    xy_ind = Floor(xcen - x_min)/covariant_bin + lookup_dim*Floor(ycen - y_min)/covariant_bin
+    ind_hist = histogram(xy_ind, min=0, max=Max(xy_ind) + 1 + lookup_dim, binsize=1, reverse_indices=ri)
+    n_offsets = 5
+    x_offsets = [0, -1, 0, 1, 0]
+    y_offsets = [-1, 0, 0, 0, 1]
+    ind_offsets = x_offsets + y_offsets*lookup_dim
+
+;    xmin = Floor(Min(xcen))
+;    ymin = Floor(Min(ycen))
+;    x_hist = histogram(xcen, binsize=covariant_bin+1, min=xmin, reverse_indices=rix)
+;    y_hist = histogram(ycen, binsize=covariant_bin+1, min=ymin, reverse_indices=riy)
+
+    FOR b_ii=0L,n_baselines-1 DO BEGIN
+        b_i = baseline_inds[b_ii]
+        bi_covariant_i = ri[ri[xy_ind[b_ii]]:ri[xy_ind[b_ii]+1]-1]
+        n_covariant = Total(ind_hist[xy_ind[b_ii] + ind_offsets])
+        n = 0
+        bi_covariant_i = Lonarr(n_covariant)
+        FOR offset_i=0,n_offsets-1 DO BEGIN
+            off_i = xy_ind[b_ii] + ind_offsets[offset_i]
+            off_n = ind_hist[off_i]
+            IF off_n GT 0 THEN BEGIN
+                bi_covariant_i[n] = ri[ri[off_i]:ri[off_i+1]-1]
+                n += off_n
+            ENDIF
+        ENDFOR
+;        ; Guaranteed at least one match, since the baseline being matched is itself in the set
+;        bi_covariant_i = where(covariance_flag, n_covariant)
+        covariances = Complexarr(n_covariant)
+
+        covariance_box = Make_array(psf_dim3, type=arr_type)
+        ; Use the zero-offset beam model as the reference
+        IF conj_flag[b_ii] THEN $
+            covariance_box[psf_base_inds] = Conj(*(*beam_arr[pol_i,fbin_use,b_i mod obs.nbaselines])[0,0]) ELSE $
+            covariance_box[psf_base_inds] = *(*beam_arr[pol_i,fbin_use,b_i mod obs.nbaselines])[0,0]
+        dx = xcen[bi_covariant_i] - xcen[b_ii]
+        dy = ycen[bi_covariant_i] - ycen[b_ii]
+        x0 = Floor(dx)
+        y0 = Floor(dy)
+        ; Use the same indexing notation as visibility_grid.pro
+        ind_offset = x0 + y0*psf_dim2
+        x_offset=Fix(Floor((dx-x0)*psf_resolution) mod psf_resolution, type=12) ; type=12 is unsigned int
+        y_offset=Fix(Floor((dy-y0)*psf_resolution) mod psf_resolution, type=12)
+
+        covariance_box2 =  Make_array(psf_dim3, n_covariant, type=arr_type)
+
+        FOR covariant_i=0L,n_covariant-1 DO BEGIN
+            b_i2 = bi_covariant_i[covariant_i]
+            beam_use = *(*beam_arr[pol_i,fbin_use,baseline_inds[b_i2] mod obs.nbaselines])[x_offset[covariant_i], y_offset[covariant_i]]
+            IF ~conj_flag[b_i2] THEN beam_use = Conj(beam_use)
+            ; use the complex conjugate for the second half of the baselines, that have been mirrored
+            ; Note that I have combined this with the below complex conjugate, so only apply Conj() to the first half
+            ; Note: this should be the complex conjugate, but I combined that and the previous complex conjugate to save computations
+            covariance_box2[covariant_i*psf_dim3 + psf_base_inds + ind_offset[covariant_i]] = beam_use
+
+            ; covariances[covariant_i] = Abs(Total(covariance_box*covariance_box2))^2./(Total(Abs(covariance_box)^2.)*Total(Abs(covariance_box2)^2.))
+            ;covariances[covariant_i] = 2*Total(covariance_box*covariance_box2)/(Total(Abs(covariance_box)^2.) + Total(Abs(covariance_box2)^2.))
+
+        ENDFOR
+        covariances=matrix_multiply((covariance_box2)/n_covariant,(covariance_box),/atranspose)
+        ija[b_ii] = Ptr_new(bi_covariant_i)
+        sa[b_ii] = Ptr_new(covariances)
+    ENDFOR
+
+    ; Save the covariances and baseline indices as a sparse matrix, in the same format as the holographic mapping function
+    covariance_map_fn = {ija:ija,sa:sa,i_use:baseline_inds,norm:obs.nbaselines,indexed:0}
+    RETURN, covariance_map_fn
+END

fhd_core/calibration/calculate_redundant_cal_correction.pro

@@ -0,0 +1,15 @@
+PRO calculate_redundant_cal_correction, vis_data, vis_model, covariance_map_fn, A_ind, B_ind, gain, freq_i=freq_i,$
+    redundant_cal_correction=redundant_cal_correction, baseline_inds=baseline_inds, covariance_threshold=covariance_threshold,$
+    redundant_calibration_weight=redundant_calibration_weight
+
+    redundant_cal_correction_use = Reform(redundant_cal_correction[freq_i, *])
+    n_baselines = N_Elements(redundant_cal_correction_use)
+    redundant_cal_delta = Complexarr(n_baselines)
+    vis_residual = vis_data*weight_invert(gain[A_ind]*Conj(gain[B_ind])) - vis_model
+    vis_delta2 = [redundant_cal_correction_use[baseline_inds], Conj(redundant_cal_correction_use[baseline_inds])]
+    vis_residual -= vis_delta2
+    IF Keyword_Set(redundant_calibration_weight) THEN vis_residual /= redundant_calibration_weight
+    SPRSAX2,covariance_map_fn,vis_residual,redundant_cal_delta,/complex
+
+    redundant_cal_correction[freq_i, *] += redundant_cal_delta/2.
+ END

fhd_core/calibration/fhd_struct_init_cal.pro

@@ -8,7 +8,9 @@ FUNCTION fhd_struct_init_cal, obs, params, skymodel, gain_arr_ptr=gain_arr_ptr,
     calibration_polyfit=calibration_polyfit, bandpass_calibrate=bandpass_calibrate, $
     cal_mode_fit=cal_mode_fit, file_path_fhd=file_path_fhd, transfer_calibration=transfer_calibration, $
     calibration_auto_initialize=calibration_auto_initialize, cal_gain_init=cal_gain_init, $
-    phase_fit_iter=phase_fit_iter, cal_amp_degree_fit=cal_amp_degree_fit,cal_phase_degree_fit=cal_phase_degree_fit,_Extra=extra
+    phase_fit_iter=phase_fit_iter, cal_amp_degree_fit=cal_amp_degree_fit,cal_phase_degree_fit=cal_phase_degree_fit,$
+    use_redundant_calibration=use_redundant_calibration, redundant_calibration_iter=redundant_calibration_iter,$
+    redundant_calibration_weight=redundant_calibration_weight, cal_covariance_threshold=cal_covariance_threshold, _Extra=extra
 
 IF N_Elements(obs) EQ 0 THEN fhd_save_io,0,obs,var='obs',/restore,file_path_fhd=file_path_fhd
 IF N_Elements(params) EQ 0 THEN fhd_save_io,0,params,var='params',/restore,file_path_fhd=file_path_fhd
@@ -38,6 +40,18 @@ ref_antenna_name=(*obs.baseline_info).tile_names[ref_antenna]
 IF N_Elements(cal_convergence_threshold) EQ 0 THEN cal_convergence_threshold=1E-7
 IF N_Elements(calibration_origin) EQ 0 THEN $
     IF Tag_exist(obs,'obsname') THEN calibration_origin=obs.obsname ELSE calibration_origin=''
+IF N_Elements(use_redundant_calibration) EQ 0 THEN use_redundant_calibration=0
+; Time averaging is not yet compatible with the redundant calibration correction
+IF Keyword_Set(use_redundant_calibration) THEN cal_time_average=0
+IF N_Elements(redundant_calibration_iter) EQ 0 THEN $
+    redundant_calibration_iter=phase_fit_iter+2 ELSE $
+    redundant_calibration_iter=Long(redundant_calibration_iter)
+IF N_Elements(cal_covariance_threshold) EQ 0 THEN $
+    cal_covariance_threshold=1E-3 ELSE $
+    cal_covariance_threshold=Float(cal_covariance_threshold)
+IF N_Elements(redundant_calibration_weight) EQ 0 THEN $
+    redundant_calibration_weight=1. ELSE $
+    redundant_calibration_weight=Float(0.<redundant_calibration_weight<1.)
 IF N_Elements(cal_gain_init) EQ 0 THEN cal_gain_init=1.
 IF N_Elements(gain_arr_ptr) EQ 0 THEN BEGIN
     gain_arr=Complexarr(n_freq,n_tile)+cal_gain_init
@@ -60,7 +74,10 @@ auto_params=Ptrarr(2)
 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,$
+    auto_initialize:auto_initialize, time_avg:cal_time_average, min_solns:min_cal_solutions,$
+    convergence_iter:[-1, -1], max_iter:max_cal_iter, phase_iter:phase_fit_iter,$
+    use_redundant:use_redundant_calibration, redundant_iter:redundant_calibration_iter,$
+    redundant_weight:redundant_calibration_weight, covariance_threshold:cal_covariance_threshold,$
     tile_A:tile_A, tile_B:tile_B, tile_names:tile_names, bin_offset:bin_offset, freq:freq, gain:gain_arr_ptr,$
     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,$

fhd_core/calibration/vis_calibrate.pro

@@ -7,7 +7,7 @@ FUNCTION vis_calibrate,vis_ptr,cal,obs,status_str,psf,params,jones,vis_weight_pt
     flag_calibration=flag_calibration,vis_model_arr=vis_model_arr,$
     calibration_auto_fit=calibration_auto_fit,cal_stop=cal_stop, model_transfer=model_transfer,$
     sim_over_calibrate=sim_over_calibrate,debug_phase_longrun=debug_phase_longrun,sim_perf_calibrate=sim_perf_calibrate,$
-    debug_ave_ref=debug_ave_ref,debug_amp_longrun=debug_amp_longrun,no_png=no_png,_Extra=extra
+    debug_ave_ref=debug_ave_ref,debug_amp_longrun=debug_amp_longrun,redundant_cal_correction=redundant_cal_correction,_Extra=extra
 
   t0_0=Systime(1)
   error=0
@@ -175,8 +175,9 @@ FUNCTION vis_calibrate,vis_ptr,cal,obs,status_str,psf,params,jones,vis_weight_pt
   FOR iter=0,calibration_flag_iterate DO BEGIN
     t2_a=Systime(1)
     IF iter LT calibration_flag_iterate THEN preserve_flag=1 ELSE preserve_flag=preserve_visibilities
-    cal=vis_calibrate_subroutine(vis_ptr,vis_model_arr,vis_weight_ptr,obs,cal,$
-      preserve_visibilities=preserve_flag,_Extra=extra)
+    redundant_cal_correction=Ptrarr(n_pol)
+    cal=vis_calibrate_subroutine(vis_ptr,vis_model_arr,vis_weight_ptr,obs,cal,psf,$
+      redundant_cal_correction=redundant_cal_correction,preserve_visibilities=preserve_flag,_Extra=extra)
     IF Keyword_Set(flag_calibration) THEN vis_calibration_flag,obs,cal,n_tile_cut=n_tile_cut,_Extra=extra
     IF Keyword_Set(n_tile_cut) THEN BREAK
   ENDFOR