Skip to content

Allow complex and double precision HEALPix #235

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
nicholebarry wants to merge 4 commits into
base: master
Choose a base branch
from
Draft

Allow complex and double precision HEALPix #235

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
6a931f3
Allow complex and double precision HEALPix
nicholebarry Sep 25, 2020
bb8153c
merge old branch
nicholebarry Oct 1, 2025
1f7e90a
make jones healpix cube
nicholebarry Oct 2, 2025
a213441
Merge branch 'master' into cross_pol_healpix
nicholebarry Oct 2, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 10 additions & 2 deletions fhd_core/HEALPix/healpix_cnv_apply.pro
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,9 +4,17 @@ t0=Systime(1)
dimension=(size(image,/dimension))[0]
elements=(size(image,/dimension))[1]

IF size(hpx_cnv,/type) EQ 10 THEN hpx_cnv_inds=(*hpx_cnv).inds ELSE hpx_cnv_inds=hpx_cnv.inds

image_vector=reform(image,Float(dimension)*elements)
IF size(hpx_cnv,/type) EQ 10 THEN hpx_map=Fltarr(N_Elements((*hpx_cnv).inds)) $
ELSE hpx_map=Fltarr(N_Elements(hpx_cnv.inds))
case size(image_vector,/type) of
4: hpx_map=fltarr(N_Elements(hpx_cnv_inds))
5: hpx_map=dblarr(N_Elements(hpx_cnv_inds))
6: hpx_map=complex(fltarr(N_Elements(hpx_cnv_inds)))
9: hpx_map=dcomplex(dblarr(N_Elements(hpx_cnv_inds)))
else: hpx_map=fltarr(N_Elements(hpx_cnv_inds))
endcase

SPRSAX2,hpx_cnv,image_vector,hpx_map,transpose=1,mask=0

timing=Systime(1)-t0
Expand Down
62 changes: 49 additions & 13 deletions fhd_core/HEALPix/healpix_snapshot_cube_generate.pro
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
PRO healpix_snapshot_cube_generate,obs_in,status_str,psf_in,cal,params,vis_arr,vis_model_arr=vis_model_arr,$
PRO healpix_snapshot_cube_generate,obs_in,status_str,psf_in,cal,params,vis_arr,vis_model_arr=vis_model_arr,jones,$
file_path_fhd=file_path_fhd,ps_dimension=ps_dimension,ps_fov=ps_fov,ps_degpix=ps_degpix,$
ps_kbinsize=ps_kbinsize,ps_kspan=ps_kspan,ps_beam_threshold=ps_beam_threshold,ps_nfreq_avg=ps_nfreq_avg,$
rephase_weights=rephase_weights,n_avg=n_avg,vis_weights=vis_weights,split_ps_export=split_ps_export,$
Expand Down Expand Up @@ -141,6 +141,14 @@ PRO healpix_snapshot_cube_generate,obs_in,status_str,psf_in,cal,params,vis_arr,v
save_uvf=save_uvf, uvf_name=uvf_name[iter],bi_use=*bi_use[iter], _Extra=extra)

t_split+=t_split1

case size(residual_arr1,/type) of
4: init_arr=fltarr(n_hpx,n_freq_use)
5: init_arr=dblarr(n_hpx,n_freq_use)
6: init_arr=complex(fltarr(n_hpx,n_freq_use))
9: init_arr=dcomplex(dblarr(n_hpx,n_freq_use))
endcase

IF dirty_flag THEN BEGIN
dirty_arr1=residual_arr1
residual_flag=0
Expand All @@ -151,15 +159,25 @@ PRO healpix_snapshot_cube_generate,obs_in,status_str,psf_in,cal,params,vis_arr,v

t_hpx0=Systime(1)

beam_squared_cube=fltarr(n_hpx,n_freq_use)
weights_cube=fltarr(n_hpx,n_freq_use)
variance_cube=fltarr(n_hpx,n_freq_use)
IF residual_flag THEN res_cube=fltarr(n_hpx,n_freq_use)
IF dirty_flag THEN dirty_cube=fltarr(n_hpx,n_freq_use)
IF model_flag THEN model_cube=fltarr(n_hpx,n_freq_use)
beam_squared_cube=real_part(init_arr)
weights_cube=init_arr
variance_cube=real_part(init_arr)
IF residual_flag THEN res_cube=init_arr
IF dirty_flag THEN dirty_cube=init_arr
IF model_flag THEN model_cube=init_arr

FOR pol_i=0,n_pol-1 DO BEGIN

IF (obs.pol_names[pol_i] EQ 'YX') then begin
If (size(init_arr,/type) EQ 6) OR (size(init_arr,/type) EQ 9) then begin
print, "Skipping YX HEALPix cube generation because it is the complex conjugate of XY"
;If complex 4-pol HEALPix are being generated, then also create the Jones HEALPix cube after the loop
; so that Stokes cubes can be generated during integration.
jones_hpx_flag = 1
continue
ENDIF
ENDIF

FOR freq_i=Long64(0),n_freq_use-1 DO BEGIN

beam_squared_cube[n_hpx*freq_i]=healpix_cnv_apply((*beam_arr[pol_i,freq_i])*nf_vis_use[freq_i],hpx_cnv)
Expand Down Expand Up @@ -187,7 +205,7 @@ PRO healpix_snapshot_cube_generate,obs_in,status_str,psf_in,cal,params,vis_arr,v
fhd_save_io,status_str,file_path_fhd=file_path_fhd,var=cube_name[iter],pol_i=pol_i,/force,_Extra=extra

ENDFOR

IF Keyword_Set(save_imagecube) THEN BEGIN
save, filename = imagecube_filepath[iter], dirty_arr1, residual_arr1, model_arr1, weights_arr1, variance_arr1, beam_arr, nf_vis_use, obs_out, /compress
ENDIF
Expand All @@ -196,10 +214,28 @@ PRO healpix_snapshot_cube_generate,obs_in,status_str,psf_in,cal,params,vis_arr,v
dirty_cube=(model_cube=(res_cube=(weights_cube=(variance_cube=(beam_squared_cube=0)))))
IF iter EQ n_iter-1 THEN undefine_fhd,beam_arr

ENDFOR
obs_out=obs ;for return
Ptr_free,vis_weights_use
timing=Systime(1)-t0
IF ~Keyword_Set(silent) THEN print,'HEALPix cube export timing: ',timing,t_split,t_hpx
ENDFOR

if keyword_set(jones_hpx_flag) then begin
; Generate Jones matrix HEALPix cube to be able to translate
; instrumental pol cubes to Stokes pol cubes during integration

jones_hpx_arr=Ptrarr(4,4,/allocate)
p_corr = Complex(FLTARR(jones.dimension,jones.elements))

FOR instr_pol=0, n_pol-1 DO BEGIN
FOR sky_pol=0, n_pol-1 DO BEGIN
p_corr[jones.inds] = *jones.jinv[instr_pol,sky_pol]
*jones_hpx_arr[instr_pol,sky_pol] = healpix_cnv_apply(p_corr,hpx_cnv)
ENDFOR
ENDFOR
save, jones_hpx_arr, filename = file_basename(file_path_fhd)+'/Healpix/'+obs.obsname+'_jones_cube.sav',jones_hpx_arr

endif

obs_out=obs ;for return
Ptr_free,vis_weights_use
timing=Systime(1)-t0
IF ~Keyword_Set(silent) THEN print,'HEALPix cube export timing: ',timing,t_split,t_hpx

END
2 changes: 1 addition & 1 deletion fhd_core/fhd_main.pro
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -293,7 +293,7 @@ IF Keyword_Set(export_images) THEN BEGIN
ENDIF

;optionally export frequency-splt Healpix cubes
IF Keyword_Set(snapshot_healpix_export) THEN healpix_snapshot_cube_generate,obs,status_str,psf,cal,params,vis_arr,$
IF Keyword_Set(snapshot_healpix_export) THEN healpix_snapshot_cube_generate,obs,status_str,psf,cal,params,vis_arr,jones,$
vis_model_arr=vis_model_arr,file_path_fhd=file_path_fhd,vis_weights=vis_weights,cmd_args=cmd_args,_Extra=extra

;Optionally fill the fhd table on the mwa_qc database located on eor-00 under the mwa username. See the python script
Expand Down
23 changes: 15 additions & 8 deletions fhd_core/gridding/vis_model_freq_split.pro
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,8 @@ FUNCTION vis_model_freq_split,obs,status_str,psf,params,vis_weights,model_uv_arr
dimension=obs.dimension
degpix=obs.degpix
residual_flag=obs.residual
;Include complex components to the image if including cross-pols
if n_pol GT 2 then complex_flag=1 else complex_flag=0

IF Min(Ptr_valid(vis_data_arr)) EQ 0 THEN BEGIN
vis_data_arr=Ptrarr(n_pol)
Expand Down Expand Up @@ -137,18 +139,22 @@ FUNCTION vis_model_freq_split,obs,status_str,psf,params,vis_weights,model_uv_arr

IF Keyword_Set(fft) THEN BEGIN
IF N_Elements(x_range)<N_Elements(y_range) GT 0 THEN BEGIN
*dirty_arr[pol_i,fi]=extract_subarray(dirty_image_generate(dirty_uv,degpix=degpix)*n_vis,x_range,y_range)
*weights_arr[pol_i,fi]=extract_subarray(dirty_image_generate(weights_holo*rephase_use,degpix=degpix)*n_vis,x_range,y_range)
*variance_arr[pol_i,fi]=extract_subarray(dirty_image_generate(variance_holo*rephase_use,degpix=degpix)*n_vis,x_range,y_range)
*dirty_arr[pol_i,fi]=$
extract_subarray(dirty_image_generate(dirty_uv,degpix=degpix,no_real=complex_flag)*n_vis,x_range,y_range)
*weights_arr[pol_i,fi]=$
extract_subarray(dirty_image_generate(weights_holo*rephase_use,degpix=degpix,no_real=complex_flag)*n_vis,x_range,y_range)
*variance_arr[pol_i,fi]=$
extract_subarray(dirty_image_generate(variance_holo*rephase_use,degpix=degpix,no_real=complex_flag)*n_vis,x_range,y_range)
IF N_Elements(model_return) GT 1 THEN $
*model_arr[pol_i,fi]=extract_subarray(dirty_image_generate(model_return,degpix=degpix)*n_vis,x_range,y_range) $
*model_arr[pol_i,fi]=$
extract_subarray(dirty_image_generate(model_return,degpix=degpix,no_real=complex_flag)*n_vis,x_range,y_range) $
ELSE IF Keyword_Set(model_flag) THEN *model_arr[pol_i,fi]=init_arr
ENDIF ELSE BEGIN
*dirty_arr[pol_i,fi]=dirty_image_generate(dirty_uv,degpix=degpix)*n_vis
*weights_arr[pol_i,fi]=dirty_image_generate(weights_holo*rephase_use,degpix=degpix)*n_vis
*variance_arr[pol_i,fi]=dirty_image_generate(variance_holo*rephase_use,degpix=degpix)*n_vis
*dirty_arr[pol_i,fi]=dirty_image_generate(dirty_uv,degpix=degpix,no_real=complex_flag)*n_vis
*weights_arr[pol_i,fi]=dirty_image_generate(weights_holo*rephase_use,degpix=degpix,no_real=complex_flag)*n_vis
*variance_arr[pol_i,fi]=dirty_image_generate(variance_holo*rephase_use,degpix=degpix,no_real=complex_flag)*n_vis
IF N_Elements(model_return) GT 1 THEN $
*model_arr[pol_i,fi]=dirty_image_generate(model_return,degpix=degpix)*n_vis $
*model_arr[pol_i,fi]=dirty_image_generate(model_return,degpix=degpix,no_real=complex_flag)*n_vis $
ELSE IF Keyword_Set(model_flag) THEN *model_arr[pol_i,fi]=init_arr
ENDELSE
ENDIF ELSE BEGIN
Expand All @@ -159,6 +165,7 @@ FUNCTION vis_model_freq_split,obs,status_str,psf,params,vis_weights,model_uv_arr
ENDELSE
IF Keyword_Set(t_grid0) THEN t_grid+=t_grid0
ENDFOR

IF ~Keyword_Set(preserve_visibilities) THEN ptr_free,vis_ptr,model_ptr

obs_out.n_vis=n_vis_use
Expand Down