diff --git a/fhd_core/calibration/apply_redundant_cal_correction.pro b/fhd_core/calibration/apply_redundant_cal_correction.pro new file mode 100644 index 000000000..10eec4831 --- /dev/null +++ b/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 diff --git a/fhd_core/calibration/calculate_baseline_covariance.pro b/fhd_core/calibration/calculate_baseline_covariance.pro new file mode 100644 index 000000000..bf0d916b7 --- /dev/null +++ b/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 diff --git a/fhd_core/calibration/calculate_redundant_cal_correction.pro b/fhd_core/calibration/calculate_redundant_cal_correction.pro new file mode 100644 index 000000000..b068e78db --- /dev/null +++ b/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 diff --git a/fhd_core/calibration/fhd_struct_init_cal.pro b/fhd_core/calibration/fhd_struct_init_cal.pro index 5a528adfd..b070d8dac 100644 --- a/fhd_core/calibration/fhd_struct_init_cal.pro +++ b/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.1 DO BEGIN vis_use=vis_data2 - vis_model_matrix=vis_model2*Conj(gain_curr[B_ind]) + IF Keyword_Set(use_redundant_calibration) THEN BEGIN + vis_model_red = apply_redundant_cal_correction(vis_model2, redundant_cal_correction=*redundant_cal_correction[pol_i], freq_i=fi,$ + baseline_inds=baseline_use[b_i_use1], /use_conjugate,$ + redundant_calibration_weight=cal.redundant_weight) + vis_model_matrix=vis_model_red*Conj(gain_curr[B_ind]) + ENDIF ELSE BEGIN + vis_model_matrix=vis_model2*Conj(gain_curr[B_ind]) + ENDELSE IF Keyword_Set(calibration_weights) THEN BEGIN FOR tile_i=0L,n_tile_use-1 DO IF n_arr[tile_i] GE min_cal_solutions THEN BEGIN xmat=vis_model_matrix[*A_ind_arr[tile_i]] @@ -180,8 +199,6 @@ FUNCTION vis_calibrate_subroutine,vis_ptr,vis_model_ptr,vis_weight_ptr,obs,cal,p BREAK ENDIF IF phase_fit_iter-i GT 0 THEN gain_new*=Abs(gain_old)*weight_invert(Abs(gain_new)) ;fit only phase at first -; IF (2.*phase_fit_iter-i GT 0) AND (phase_fit_iter-i LE 0) THEN $ -; gain_new*=Mean(Abs(gain_new[where(gain_new)]))*weight_invert(Abs(gain_new)) ;then fit only average amplitude gain_curr=(gain_new+gain_old)/2. dgain=Abs(gain_curr)*weight_invert(Abs(gain_old)) diverge_i=where(dgain LT Abs(gain_old)/2.,n_diverge) @@ -190,6 +207,11 @@ FUNCTION vis_calibrate_subroutine,vis_ptr,vis_model_ptr,vis_weight_ptr,obs,cal,p if ~keyword_set(no_ref_tile) then begin gain_curr*=Conj(gain_curr[ref_tile_use])/Abs(gain_curr[ref_tile_use]) endif + IF Keyword_Set(use_redundant_calibration) AND (redundant_fit_iter-i GT 0) THEN BEGIN + calculate_redundant_cal_correction,vis_data2, vis_model2, covariance_map_fn, A_ind, B_ind, gain_old, freq_i=fi,$ + redundant_cal_correction=*redundant_cal_correction[pol_i], baseline_inds=baseline_use[b_i_use1],$ + covariance_threshold=cal.covariance_threshold, redundant_calibration_weight=cal.redundant_weight + ENDIF conv_test[fii,i]=Max(Abs(gain_curr-gain_old)*weight_invert(Abs(gain_old))) IF i GT phase_fit_iter THEN IF conv_test[fii,i] LE conv_thresh THEN BEGIN n_converged += 1 @@ -216,6 +238,7 @@ FUNCTION vis_calibrate_subroutine,vis_ptr,vis_model_ptr,vis_weight_ptr,obs,cal,p *cal_return.gain[pol_i]=gain_arr cal_return.convergence[pol_i]=Ptr_new(convergence) cal_return.n_converged[pol_i] = n_converged + IF Tag_exist(cal_return, "convergence_iter") THEN cal_return.convergence_iter[pol_i] = i ENDFOR vis_count_i=where(weight,n_vis_cal) diff --git a/fhd_core/fhd_main.pro b/fhd_core/fhd_main.pro index 53bcf1cea..276e86b4d 100644 --- a/fhd_core/fhd_main.pro +++ b/fhd_core/fhd_main.pro @@ -44,6 +44,7 @@ IF data_flag LE 0 THEN BEGIN uvfits_read,hdr,params,layout,vis_arr,vis_weights,file_path_vis=file_path_vis,n_pol=n_pol,silent=silent,error=error,_Extra=extra IF Keyword_Set(error) THEN BEGIN print,"Error occured while reading uvfits data. Returning." + IF Keyword_Set(!Journal) THEN Journal ;write and close log file if present RETURN ENDIF fhd_save_io,status_str,layout,var='layout',/compress,file_path_fhd=file_path_fhd,_Extra=extra ;save layout structure right away for debugging. Will be overwritten a few times before the end @@ -91,6 +92,7 @@ IF data_flag LE 0 THEN BEGIN flag_calibration=flag_calibration,_Extra=extra IF Keyword_Set(error) THEN BEGIN print,"Error occured while attempting to transfer weights. Returning." + IF Keyword_Set(!Journal) THEN Journal ;write and close log file if present RETURN ENDIF @@ -121,6 +123,7 @@ IF data_flag LE 0 THEN BEGIN ; This is checked after optionally transfering or modifying the weights, and after writing the settings file. error=1 print,"All tiles flagged! No data left to use. Returning" + IF Keyword_Set(!Journal) THEN Journal ;write and close log file if present RETURN ENDIF @@ -142,12 +145,14 @@ IF data_flag LE 0 THEN BEGIN transfer_calibration=transfer_calibration,timing=cal_timing,error=error,model_uv_arr=model_uv_arr,$ return_cal_visibilities=return_cal_visibilities,vis_model_arr=vis_model_arr,$ calibration_visibilities_subtract=calibration_visibilities_subtract,silent=silent,$ - flag_calibration=flag_calibration,cal_stop=cal_stop,_Extra=extra) + flag_calibration=flag_calibration,cal_stop=cal_stop, + redundant_cal_correction=redundant_cal_correction,_Extra=extra) IF Keyword_Set(return_cal_visibilities) THEN $ vis_calibrate_qu_mixing,vis_arr,vis_model_arr,vis_weights,obs,cal IF ~Keyword_Set(silent) THEN print,String(format='("Calibration timing: ",A)',Strn(cal_timing)) IF Keyword_Set(error) THEN BEGIN print,"Error occured during calibration. Returning." + IF Keyword_Set(!Journal) THEN Journal ;write and close log file if present RETURN ENDIF fhd_save_io,status_str,cal,var='cal',/compress,file_path_fhd=file_path_fhd,_Extra=extra @@ -158,6 +163,7 @@ IF data_flag LE 0 THEN BEGIN run_report, start_mem, t0, silent=silent RETURN endif + redundant_correction_flag = cal.use_redundant ENDIF IF Keyword_Set(flag_visibilities) THEN BEGIN @@ -191,6 +197,7 @@ IF data_flag LE 0 THEN BEGIN IF Keyword_Set(error) THEN BEGIN print,"Error occured during source modeling" + IF Keyword_Set(!Journal) THEN Journal ;write and close log file if present RETURN ENDIF @@ -220,6 +227,7 @@ IF data_flag LE 0 THEN BEGIN t_save0=Systime(1) vis_export,obs,status_str,vis_arr,vis_weights,file_path_fhd=file_path_fhd,/compress IF Keyword_Set(model_flag) THEN vis_export,obs,status_str,vis_model_arr,vis_weights,file_path_fhd=file_path_fhd,/compress,/model + IF Keyword_Set(redundant_correction_flag) THEN vis_export,obs,status_str,redundant_cal_correction,vis_weights,file_path_fhd=file_path_fhd,/compress,/redundant_correction t_save=Systime(1)-t_save0 IF ~Keyword_Set(silent) THEN print,'Visibility save time: ',t_save ENDIF @@ -228,8 +236,10 @@ IF data_flag LE 0 THEN BEGIN IF Keyword_Set(grid_recalculate) THEN BEGIN image_uv_arr=visibility_grid_wrap(vis_arr,vis_weights,obs,status_str,psf,params,file_path_fhd=file_path_fhd,vis_model_arr=vis_model_arr,$ deconvolve=deconvolve,model_flag=model_flag,snapshot_healpix_export=snapshot_healpix_export,mapfn_recalculate=mapfn_recalculate,$ - save_visibilities=save_visibilities,error=error,no_save=no_save,weights_arr=weights_arr,model_uv_holo=model_uv_holo,$ - return_decon_visibilities=return_decon_visibilities,_Extra=extra) + save_visibilities=save_visibilities,error=error,no_save=no_save,weights_arr=weights_arr,model_uv_arr=model_uv_arr,$ + return_decon_visibilities=return_decon_visibilities,$ + redundant_correction_flag=redundant_correction_flag,redundant_cal_correction=redundant_cal_correction,redundantCorr_uv_arr=redundantCorr_uv_arr,$ + _Extra=extra) ENDIF ELSE BEGIN print,'Visibilities not re-gridded' ENDELSE @@ -255,7 +265,7 @@ ENDELSE ;Generate fits data files and images IF Keyword_Set(export_images) THEN BEGIN fhd_quickview,obs,status_str,psf,cal,jones,skymodel,fhd_params,image_uv_arr=image_uv_arr,weights_arr=weights_arr,$ - model_uv_holo=model_uv_holo,beam_arr=beam_arr,file_path_fhd=file_path_fhd,silent=silent,$ + model_uv_arr=model_uv_arr,redundantCorr_uv_arr=redundantCorr_uv_arr,beam_arr=beam_arr,file_path_fhd=file_path_fhd,silent=silent,$ map_fn_arr=map_fn_arr,transfer_mapfn=transfer_mapfn,_Extra=extra ENDIF diff --git a/fhd_core/gridding/visibility_grid.pro b/fhd_core/gridding/visibility_grid.pro index a65c427ed..c1aa35912 100644 --- a/fhd_core/gridding/visibility_grid.pro +++ b/fhd_core/gridding/visibility_grid.pro @@ -6,6 +6,7 @@ FUNCTION visibility_grid,visibility_ptr,vis_weight_ptr,obs,status_str,psf,params model_ptr=model_ptr,model_return=model_return,preserve_visibilities=preserve_visibilities,$ error=error,grid_uniform=grid_uniform,$ grid_spectral=grid_spectral,spectral_uv=spectral_uv,spectral_model_uv=spectral_model_uv,$ + redundantCorr_return=redundantCorr_return, redundantCorr_ptr=redundantCorr_ptr,$ beam_per_baseline=beam_per_baseline,uv_grid_phase_only=uv_grid_phase_only,_Extra=extra t0_0=Systime(1) heap_gc @@ -73,6 +74,20 @@ IF Ptr_valid(model_ptr) THEN BEGIN ; ELSE vis_arr_use-=(Temporary(*model_ptr))[vis_inds_use] ENDELSE ENDIF +redundantCorr_flag=0 +IF Ptr_valid(redundantCorr_ptr) THEN BEGIN + IF Keyword_Set(redundantCorr_return) THEN BEGIN + IF Keyword_Set(preserve_visibilities) THEN redundantCorr_use=(*redundantCorr_ptr)[vis_inds_use] ELSE BEGIN + redundantCorr_use=(Temporary(*redundantCorr_ptr))[vis_inds_use] + ptr_free,redundantCorr_ptr + ENDELSE + IF Keyword_Set(double_precision) THEN redundantCorr_return=DComplexarr(dimension,elements) ELSE redundantCorr_return=Complexarr(dimension,elements) + redundantCorr_flag=1 + ENDIF ELSE BEGIN +; IF Keyword_Set(preserve_visibilities) THEN vis_arr_use-=(*model_ptr)[vis_inds_use] $ +; ELSE vis_arr_use-=(Temporary(*model_ptr))[vis_inds_use] + ENDELSE +ENDIF frequency_array=(*obs.baseline_info).freq frequency_array=frequency_array[fi_use] @@ -140,6 +155,7 @@ IF n_conj GT 0 THEN BEGIN ww[conj_i]=-ww[conj_i] vis_arr_use[*,conj_i]=Conj(vis_arr_use[*,conj_i]) IF model_flag THEN model_use[*,conj_i]=Conj(model_use[*,conj_i]) + IF redundantCorr_flag THEN redundantCorr_use[*,conj_i]=Conj(redundantCorr_use[*,conj_i]) ENDIF xcen=Float(frequency_array#Temporary(kx_arr)) ycen=Float(frequency_array#Temporary(ky_arr)) @@ -316,6 +332,10 @@ FOR bi=0L,n_bin_use-1 DO BEGIN model_box1=model_use[inds] model_box=model_box1[xyf_ui] ENDIF + IF redundantCorr_flag THEN BEGIN + redundantCorr_box1=redundantCorr_use[inds] + redundantCorr_box=redundantCorr_box1[xyf_ui] + ENDIF repeat_i=where(psf_weight GT 1,n_rep,complement=single_i,ncom=n_single) @@ -326,11 +346,14 @@ FOR bi=0L,n_bin_use-1 DO BEGIN IF model_flag THEN FOR rep_ii=0,n_rep-1 DO $ model_box[repeat_i[rep_ii]]=Total(model_box1[xyf_ui0[rep_ii]:xyf_ui[rep_ii]]) + IF redundantCorr_flag THEN FOR rep_ii=0,n_rep-1 DO $ + redundantCorr_box[repeat_i[rep_ii]]=Total(redundantCorr_box1[xyf_ui0[rep_ii]:xyf_ui[rep_ii]]) vis_n=n_xyf_bin ENDIF ELSE BEGIN rep_flag=0 IF model_flag THEN model_box=model_use[inds] + IF redundantCorr_flag THEN redundantCorr_box=redundantCorr_use[inds] vis_box=vis_arr_use[inds] psf_weight=Replicate(1.,vis_n) bt_index = inds / n_freq_use @@ -391,6 +414,10 @@ FOR bi=0L,n_bin_use-1 DO BEGIN box_arr=matrix_multiply(Temporary(model_box)/n_vis,box_matrix_dag,/atranspose,/btranspose) model_return[xmin_use:xmin_use+psf_dim-1,ymin_use:ymin_use+psf_dim-1]+=Temporary(box_arr) ENDIF + IF redundantCorr_flag THEN BEGIN + box_arr=matrix_multiply(Temporary(redundantCorr_box)/n_vis,box_matrix_dag,/atranspose,/btranspose) + redundantCorr_return[xmin_use:xmin_use+psf_dim-1,ymin_use:ymin_use+psf_dim-1]+=Temporary(box_arr) + ENDIF box_arr=matrix_multiply(Temporary(vis_box)/n_vis,box_matrix_dag,/atranspose,/btranspose) image_uv[xmin_use:xmin_use+psf_dim-1,ymin_use:ymin_use+psf_dim-1]+=Temporary(box_arr) @@ -465,6 +492,7 @@ IF Keyword_Set(grid_uniform) THEN BEGIN IF weights_flag THEN weights*=weight_invert(filter_use) IF variance_flag THEN variance*=weight_invert(filter_use) IF model_flag THEN model_return*=weight_invert(filter_use) + IF redundantCorr_flag THEN redundantCorr_return*=weight_invert(filter_use) ENDIF IF ~Keyword_Set(no_conjugate) THEN BEGIN @@ -472,6 +500,7 @@ IF ~Keyword_Set(no_conjugate) THEN BEGIN IF weights_flag THEN weights=(weights+conjugate_mirror(weights))/2. IF variance_flag THEN variance=(variance+conjugate_mirror(variance))/4. ;2? IF model_flag THEN model_return=(model_return+conjugate_mirror(model_return))/2. + IF redundantCorr_flag THEN redundantCorr_return=(redundantCorr_return+conjugate_mirror(redundantCorr_return))/2. IF uniform_flag THEN uniform_filter=(uniform_filter+conjugate_mirror(uniform_filter))/2. ENDIF diff --git a/fhd_core/gridding/visibility_grid_wrap.pro b/fhd_core/gridding/visibility_grid_wrap.pro index cfbe4f3b8..30189b8f5 100644 --- a/fhd_core/gridding/visibility_grid_wrap.pro +++ b/fhd_core/gridding/visibility_grid_wrap.pro @@ -1,6 +1,8 @@ FUNCTION visibility_grid_wrap,vis_arr,vis_weights,obs,status_str,psf,params,file_path_fhd=file_path_fhd,vis_model_arr=vis_model_arr,$ deconvolve=deconvolve,model_flag=model_flag,snapshot_healpix_export=snapshot_healpix_export,mapfn_recalculate=mapfn_recalculate,$ - save_visibilities=save_visibilities,error=error,no_save=no_save,weights_arr=weights_arr,model_uv_holo=model_uv_holo,_Extra=extra + save_visibilities=save_visibilities,error=error,no_save=no_save,weights_arr=weights_arr,model_uv_arr=model_uv_arr,$ + redundant_correction_flag=redundant_correction_flag, redundant_cal_correction=redundant_cal_correction,redundantCorr_uv_arr=redundantCorr_uv_arr,$ + _Extra=extra n_pol=obs.n_pol t_grid=fltarr(n_pol) @@ -10,10 +12,12 @@ IF Keyword_Set(deconvolve) THEN map_fn_arr=Ptrarr(n_pol) image_uv_arr=Ptrarr(n_pol,/allocate) weights_arr=Ptrarr(n_pol,/allocate) -IF Keyword_Set(model_flag) THEN model_uv_holo=Ptrarr(n_pol,/allocate) +IF Keyword_Set(model_flag) THEN model_uv_arr=Ptrarr(n_pol,/allocate) +IF Keyword_Set(redundant_correction_flag) THEN redundantCorr_uv_arr=Ptrarr(n_pol,/allocate) IF N_Elements(weights_grid) EQ 0 THEN weights_grid=1 FOR pol_i=0,n_pol-1 DO BEGIN IF Keyword_Set(model_flag) THEN model_return=1 + IF Keyword_Set(redundant_correction_flag) THEN redundantCorr_return=1 IF Keyword_Set(snapshot_healpix_export) THEN preserve_visibilities=1 ELSE preserve_visibilities=0 IF Keyword_Set(preserve_visibilities) THEN return_mapfn=0 ELSE return_mapfn=mapfn_recalculate IF Keyword_Set(mapfn_recalculate) AND Keyword_Set(save_visibilities) THEN preserve_vis_grid=0 ELSE preserve_vis_grid=preserve_visibilities @@ -22,8 +26,9 @@ FOR pol_i=0,n_pol-1 DO BEGIN grid_uv=visibility_grid(vis_arr[pol_i],vis_weights[pol_i],obs,status_str,psf,params,file_path_fhd=file_path_fhd,$ timing=t_grid0,polarization=pol_i,weights=weights_grid,silent=silent,uniform_filter=uniform_filter,$ mapfn_recalculate=mapfn_recalculate,return_mapfn=return_mapfn,error=error,no_save=no_save,$ - model_return=model_return,model_ptr=vis_model_arr[pol_i],preserve_visibilities=preserve_vis_grid,$ - no_conjugate=no_conjugate,_Extra=extra) + model_return=model_return,model_ptr=vis_model_arr[pol_i],$ + redundantCorr_return=redundantCorr_return, redundantCorr_ptr=redundant_cal_correction[pol_i],$ + preserve_visibilities=preserve_vis_grid, no_conjugate=no_conjugate,_Extra=extra) IF Keyword_Set(error) THEN BEGIN print,"Error occured during gridding. Returning." IF Keyword_Set(!Journal) THEN Journal ;write and close log file if present @@ -35,7 +40,9 @@ FOR pol_i=0,n_pol-1 DO BEGIN IF Keyword_Set(deconvolve) THEN IF Keyword_Set(return_mapfn) THEN map_fn_arr[pol_i]=Ptr_new(return_mapfn,/no_copy) *image_uv_arr[pol_i]=Temporary(grid_uv) IF Keyword_Set(model_flag) THEN $ - *model_uv_holo[pol_i]=Temporary(model_return) + *model_uv_arr[pol_i]=Temporary(model_return) + IF Keyword_Set(redundant_correction_flag) THEN $ + *redundantCorr_uv_arr[pol_i]=Temporary(redundantCorr_return) IF N_Elements(weights_grid) GT 0 THEN BEGIN *weights_arr[pol_i]=Temporary(weights_grid) weights_grid=1 @@ -45,15 +52,17 @@ ENDFOR image_uv_arr=crosspol_reformat(image_uv_arr) weights_arr=crosspol_reformat(weights_arr) IF Keyword_Set(model_flag) THEN $ - model_uv_holo=crosspol_reformat(model_uv_holo) + model_uv_arr=crosspol_reformat(model_uv_arr) FOR pol_i=0, n_pol-1 DO BEGIN fhd_save_io,status_str,*image_uv_arr[pol_i],var='grid_uv',/compress,file_path_fhd=file_path_fhd,pol_i=pol_i,obs=obs,_Extra=extra fhd_save_io,status_str,*weights_arr[pol_i],var='weights_uv',/compress,file_path_fhd=file_path_fhd,pol_i=pol_i,obs=obs,_Extra=extra IF Keyword_Set(model_flag) THEN $ - fhd_save_io,status_str,*model_uv_holo[pol_i],var='grid_uv_model',/compress,file_path_fhd=file_path_fhd,pol_i=pol_i,obs=obs,_Extra=extra + fhd_save_io,status_str,*model_uv_arr[pol_i],var='grid_uv_model',/compress,file_path_fhd=file_path_fhd,pol_i=pol_i,obs=obs,_Extra=extra + IF Keyword_Set(redundant_correction_flag) THEN $ + fhd_save_io,status_str,*redundantCorr_uv_arr[pol_i],var='grid_uv_redundantcorr',/compress,file_path_fhd=file_path_fhd,pol_i=pol_i,obs=obs,_Extra=extra ENDFOR fhd_save_io,status_str,obs,var='obs',/compress,file_path_fhd=file_path_fhd,_Extra=extra ;need to save here to enter nf_vis IF ~Keyword_Set(silent) THEN print,'Gridding time:',t_grid RETURN,image_uv_arr -END \ No newline at end of file +END diff --git a/fhd_core/setup_metadata/fhd_save_io.pro b/fhd_core/setup_metadata/fhd_save_io.pro index 911618db9..501b36f7b 100644 --- a/fhd_core/setup_metadata/fhd_save_io.pro +++ b/fhd_core/setup_metadata/fhd_save_io.pro @@ -39,7 +39,8 @@ ENDIF IF Keyword_Set(restore) THEN no_save=1 IF Keyword_Set(reset) THEN status_str={hdr:0,params:0,obs:0,layout:0,psf:0,antenna:0,jones:0,cal:0,skymodel:0,source_array:0,vis_weights:0,auto_corr:0,$ - vis_ptr:intarr(4),vis_model_ptr:intarr(4),grid_uv:intarr(4),weights_uv:intarr(4),grid_uv_model:intarr(4),vis_count:0,$ + vis_ptr:intarr(4),vis_model_ptr:intarr(4),vis_redundantCorr_ptr:intarr(4),$ + grid_uv:intarr(4),weights_uv:intarr(4),grid_uv_model:intarr(4),grid_uv_redundantCorr:intarr(4),vis_count:0,$ map_fn:intarr(4),fhd:0,fhd_params:0,hpx_cnv:0,healpix_cube:intarr(4),hpx_even:intarr(4),hpx_odd:intarr(4),complete:0} IF size(status_str,/type) NE 8 THEN status_str=getvar_savefile(status_path+'.sav','status_str', compatibility_mode = compatibility_mode) status_use=status_str @@ -70,9 +71,11 @@ CASE var_name OF ;listed in order typically generated 'auto_corr':BEGIN status_use.auto_corr=1 & path_add='_autos' & subdir='vis_data' & obs_flag=1 & END 'vis_ptr':BEGIN status_use.vis_ptr[pol_i]=1 & path_add='_vis_'+pol_names[pol_i] & subdir='vis_data' & obs_flag=1 & pol_flag=1 & END 'vis_model_ptr':BEGIN status_use.vis_model_ptr[pol_i]=1 & path_add='_vis_model_'+pol_names[pol_i] & subdir='vis_data' & obs_flag=1 & pol_flag=1 & END + 'vis_redundantcorr_ptr':BEGIN status_use.vis_redundantCorr_ptr[pol_i]=1 & path_add='_vis_model_'+pol_names[pol_i] & subdir='vis_data' & obs_flag=1 & pol_flag=1 & END 'grid_uv':BEGIN status_use.grid_uv[pol_i]=1 & path_add='_uv_'+pol_names[pol_i] & subdir='grid_data'& pol_flag=1 & END 'weights_uv':BEGIN status_use.weights_uv[pol_i]=1 & path_add='_uv_weights_'+pol_names[pol_i] & subdir='grid_data'& pol_flag=1 & END 'grid_uv_model':BEGIN status_use.grid_uv_model[pol_i]=1 & path_add='_uv_model_'+pol_names[pol_i] & subdir='grid_data'& pol_flag=1 & END + 'grid_uv_redundantcorr':BEGIN status_use.grid_uv_redundantCorr[pol_i]=1 & path_add='_uv_grid_uv_redundantCorr_'+pol_names[pol_i] & subdir='grid_data'& pol_flag=1 & END 'vis_count':BEGIN status_use.vis_count=1 & path_add='_vis_count' & subdir='grid_data'& END ;NOTE: Because of it's size, only the map_fn can be saved in the mapfn .sav file. It has to be restored using RESTORE , ; and including other parameters can cause unwanted behavior when they are restored @@ -137,4 +140,4 @@ IF Keyword_Set(status_save) THEN BEGIN IF file_test(dir_use) EQ 0 THEN file_mkdir,dir_use SAVE,status_str,filename=status_path+'.sav' ENDIF -END \ No newline at end of file +END diff --git a/fhd_output/fhd_quickview.pro b/fhd_output/fhd_quickview.pro index 0acf4f99e..db92e07b8 100644 --- a/fhd_output/fhd_quickview.pro +++ b/fhd_output/fhd_quickview.pro @@ -1,5 +1,5 @@ PRO fhd_quickview,obs,status_str,psf,cal,jones,skymodel,fhd_params,image_uv_arr=image_uv_arr,weights_arr=weights_arr,$ - model_uv_arr=model_uv_arr,file_path_fhd=file_path_fhd,silent=silent,show_grid=show_grid,$ + model_uv_arr=model_uv_arr,redundantCorr_uv_arr=redundantCorr_uv_arr,file_path_fhd=file_path_fhd,silent=silent,show_grid=show_grid,$ gridline_image_show=gridline_image_show,pad_uv_image=pad_uv_image,image_filter_fn=image_filter_fn,$ grid_spacing=grid_spacing,reverse_image=reverse_image,show_obsname=show_obsname,mark_zenith=mark_zenith,$ no_fits=no_fits,no_png=no_png,ring_radius=ring_radius,zoom_low=zoom_low,zoom_high=zoom_high,zoom_radius=zoom_radius,$ @@ -91,9 +91,20 @@ IF N_Elements(model_uv_arr) EQ 0 THEN BEGIN model_uv_arr=getvar_savefile(fhd_sav_filepath+'.sav','model_uv_holo') ENDIF ELSE model_flag=0 ENDIF +IF N_Elements(redundantCorr_uv_arr) EQ 0 THEN BEGIN + IF Min(status_str.grid_uv_redundantCorr[0:n_pol-1]) GT 0 THEN BEGIN + redundantCorr_uv_arr=Ptrarr(n_pol,/allocate) + FOR pol_i=0,n_pol-1 DO BEGIN + fhd_save_io,status_str,grid_uv_redundantCorr,var='grid_uv_redundantCorr',/restore,file_path_fhd=file_path_fhd,obs=obs,pol_i=pol_i,_Extra=extra + *redundantCorr_uv_arr[pol_i]=grid_uv_redundantCorr + ENDFOR + ENDIF +ENDIF IF residual_flag THEN model_flag=0 +redundantCorr_flag = Max(Ptr_valid(redundantCorr_uv_arr)) + IF Keyword_Set(image_filter_fn) THEN BEGIN dummy_img=Call_function(image_filter_fn,fltarr(2,2),name=filter_name,/return_name_only) IF Keyword_Set(filter_name) THEN filter_name='_'+filter_name ELSE filter_name='' @@ -261,6 +272,7 @@ instr_dirty_arr=Ptrarr(n_pol) instr_sources=Ptrarr(n_pol) instr_rings=Ptrarr(n_pol) filter_arr=Ptrarr(n_pol,/allocate) +instr_redundantCorr_arr = Ptrarr(n_pol) FOR pol_i=0,n_pol-1 DO BEGIN complex_flag = pol_i GT 1 ? 1 : 0 instr_dirty_arr[pol_i]=Ptr_new(dirty_image_generate(*image_uv_arr[pol_i],degpix=degpix,weights=*weights_arr[pol_i],/antialias,$ @@ -269,6 +281,8 @@ FOR pol_i=0,n_pol-1 DO BEGIN IF model_flag THEN instr_model_arr[pol_i]=Ptr_new(dirty_image_generate(*model_uv_arr[pol_i],degpix=degpix,weights=*weights_arr[pol_i],/antialias,$ image_filter_fn=image_filter_fn,pad_uv_image=pad_uv_image,file_path_fhd=file_path_fhd,filter=filter_arr[pol_i],$ beam_ptr=beam_base_out[pol_i],no_real=complex_flag,_Extra=extra)) + IF redundantCorr_flag THEN instr_redundantCorr_arr[pol_i]=Ptr_new(dirty_image_generate(*redundantCorr_uv_arr[pol_i],degpix=degpix,weights=*weights_arr[pol_i],/antialias,$ + image_filter_fn=image_filter_fn,pad_uv_image=pad_uv_image,file_path_fhd=file_path_fhd,filter=filter_arr[pol_i],beam_ptr=beam_base_out[pol_i],_Extra=extra)) IF source_flag THEN BEGIN IF Keyword_Set(ring_radius) THEN instr_rings[pol_i]=Ptr_new(source_image_generate(source_arr_out,obs_out,pol_i=pol_i,resolution=16,$ dimension=dimension,restored_beam_width=restored_beam_width,ring_radius=ring_radius,_Extra=extra)) @@ -287,6 +301,7 @@ for pol_i=0,n_pol-1 do begin IF model_flag THEN *instr_model_arr[pol_i]*=renorm_factor[pol_i] IF source_flag THEN BEGIN *instr_sources[pol_i]*=renorm_factor[pol_i] + IF redundantCorr_flag THEN *instr_redundantCorr_arr[pol_i]*=renorm_factor[pol_i] IF Keyword_Set(ring_radius) THEN *instr_rings[pol_i]*=renorm_factor[pol_i] ENDIF endfor @@ -421,6 +436,7 @@ FOR pol_i=0,n_pol-1 DO BEGIN stokes_model=(*stokes_model_arr[pol_i])*beam_mask stokes_dirty=(*stokes_dirty_arr[pol_i])*beam_mask ENDIF + IF redundantCorr_flag THEN instr_redundantCorr=*instr_redundantCorr_arr[pol_i]*(*beam_correction_out[pol_i]) stokes_residual=(*stokes_residual_arr[pol_i])*beam_mask IF source_flag THEN BEGIN instr_source=*instr_sources[pol_i] @@ -476,6 +492,14 @@ FOR pol_i=0,n_pol-1 DO BEGIN offset_lat=offset_lat,offset_lon=offset_lon,label_spacing=label_spacing,map_reverse=map_reverse,show_grid=show_grid,$ title=title_fhd,/sphere,astr=astr_out2,contour_image=beam_contour_stokes,_Extra=extra ENDIF + IF redundantCorr_flag THEN BEGIN + redundantCorr_low_use=Min(instr_redundantCorr[beam_i])>(-5.*Stddev(instr_redundantCorr[beam_i])) + redundantCorr_high_use=Max(instr_redundantCorr[beam_i])<(10.*Stddev(instr_redundantCorr[beam_i])) + + Imagefast,instr_redundantCorr[zoom_low:zoom_high,zoom_low:zoom_high]+mark_image,file_path=image_path+filter_name+'_RedundantCorrection_'+pol_names[pol_i],$ + /right,sig=2,color_table=0,back='white',reverse_image=reverse_image,low=redundantCorr_low_use,high=redundantCorr_high_use,$ + title=title_fhd,show_grid=show_grid,astr=astr_out2,contour_image=beam_contour_arr[pol_i],_Extra=extra + ENDIF Imagefast,instr_residual[zoom_low:zoom_high,zoom_low:zoom_high]+mark_image,file_path=image_path+filter_name+res_name+pol_names[pol_i],$ /right,sig=2,color_table=0,back='white',reverse_image=reverse_image,low=instr_low_use,high=instr_high_use,$ title=title_fhd,show_grid=show_grid,astr=astr_out2,contour_image=beam_contour_arr[pol_i],_Extra=extra @@ -507,6 +531,9 @@ FOR pol_i=0,n_pol-1 DO BEGIN FitsFast,stokes_model,fits_header_apparent,/write,file_path=output_path+filter_name+'_Model_'+pol_names[pol_i+4] FitsFast,stokes_dirty,fits_header_apparent,/write,file_path=output_path+filter_name+'_Dirty_'+pol_names[pol_i+4] ENDIF + IF redundantCorr_flag THEN BEGIN + FitsFast,instr_redundantCorr,fits_header_apparent,/write,file_path=output_path+filter_name+'_RedundantCorrection_'+pol_names[pol_i] + ENDIF FitsFast,instr_residual,fits_header_apparent,/write,file_path=output_path+filter_name+res_name+pol_names[pol_i] FitsFast,beam_use,fits_header,/write,file_path=output_path+'_Beam_'+pol_names[pol_i] IF weights_flag THEN FitsFast,Abs(*weights_use[pol_i])*obs.n_vis,fits_header_uv,/write,file_path=output_path+'_UV_weights_'+pol_names[pol_i] diff --git a/fhd_output/vis_export.pro b/fhd_output/vis_export.pro index c4f7024c6..bfff5cf01 100644 --- a/fhd_output/vis_export.pro +++ b/fhd_output/vis_export.pro @@ -1,13 +1,19 @@ -PRO vis_export,obs,status_str,vis_ptr_arr,vis_weight_ptr,file_path_fhd=file_path_fhd,pol_i=pol_i,compress=compress,model=model +PRO vis_export,obs,status_str,vis_ptr_arr,vis_weight_ptr,file_path_fhd=file_path_fhd,pol_i=pol_i,compress=compress,$ + model=model, redundant_correction=redundant_correction IF N_Elements(compress) EQ 0 THEN compress=1 pol_names=obs.pol_names res_name='Residual' +var_name='vis_ptr' IF obs.residual EQ 0 THEN res_name='Dirty' IF Keyword_Set(model) THEN BEGIN res_name='Model' var_name='vis_model_ptr' -ENDIF ELSE var_name='vis_ptr' +ENDIF +IF Keyword_Set(redundant_correction) THEN BEGIN + res_name = 'Redundant calibration correction' + var_name = 'vis_redundantCorr_ptr' +ENDIF IF min(Ptr_valid(vis_ptr_arr)) EQ 0 THEN BEGIN print,res_name+" visibilities NULL! Not exported!"