Skip to content
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

Interpolation test #232

Draft
wants to merge 3 commits into
base: master
Choose a base branch
from
+89 −24
Draft

Interpolation test #232

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
5f6e6d8
Interpolation test commits"
lmberkhout Jul 8, 2020
5be1031
Merge branch 'master' of https://github.com/EoRImaging/FHD into inter…
lmberkhout Jul 8, 2020
24a285d
commits for review
lmberkhout Sep 18, 2020
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
48 changes: 31 additions & 17 deletions fhd_core/HEALPix/healpix_cnv_generate.pro
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ elements=obs.elements

IF N_Elements(hpx_radius) EQ 0 THEN BEGIN
IF Keyword_Set(mask) THEN BEGIN

xv_arr=meshgrid(dimension,elements,1)
yv_arr=meshgrid(dimension,elements,2)
;set /ignore_refraction for speed, since we don't need to be exact
Expand All @@ -32,15 +33,14 @@ ENDIF ELSE radius=hpx_radius

;all angles in DEGREES
;uses RING index scheme

;check if a string, if it is assume it is a filepath to a save file with the desired indices
; (will NOT be over-written with the indices)
IF Keyword_Set(restrict_hpx_inds) AND (size(restrict_hpx_inds,/type) NE 7) THEN restrict_hpx_inds=observation_healpix_inds_select(obs)
IF size(restrict_hpx_inds,/type) EQ 7 THEN BEGIN
file_path_use=restrict_hpx_inds
IF file_test(file_path_use) EQ 0 THEN file_path_use=filepath(file_path_use,root=Rootdir('fhd'),subdir='Observations')
IF file_test(file_path_use) EQ 0 THEN file_path_use=filepath(file_path_use,root=Rootdir('fhd'),subdir='Observations')

IF file_test(file_path_use) THEN BEGIN
IF file_test(file_path_use) THEN BEGIN
hpx_inds=getvar_savefile(file_path_use,'hpx_inds')
nside_test=getvar_savefile(file_path_use,names=sav_contents)
IF Max(strmatch(StrLowCase(sav_contents),'nside')) EQ 1 THEN nside=getvar_savefile(file_path_use,'nside') ELSE BEGIN
Expand All @@ -52,6 +52,7 @@ ENDIF
IF ~Keyword_Set(nside) THEN BEGIN
pix_sky=4.*!Pi*!RaDeg^2./Product(Abs(obs.astr.cdelt))
Nside=2.^(Ceil(ALOG(Sqrt(pix_sky/12.))/ALOG(2))) ;=1024. for 0.1119 degrees/pixel

ENDIF
npix=nside2npix(nside)

Expand All @@ -60,9 +61,18 @@ IF Keyword_Set(divide_pixel_area) THEN BEGIN
ENDIF ELSE pixel_area_cnv=1. ;turn this off for now

IF N_Elements(hpx_inds) GT 1 THEN BEGIN
pix2vec_ring,nside,hpx_inds,pix_coords
vec2ang,pix_coords,pix_dec,pix_ra,/astro
apply_astrometry, obs, ra_arr=pix_ra, dec_arr=pix_dec, x_arr=xv_hpx, y_arr=yv_hpx, /ad2xy
;set up orthoslant grid
x_vals = repmat(findgen(obs.dimension), 1, obs.dimension)
y_vals = Transpose(x_vals)
;shape them correctly
x_vals = reform(x_vals, obs.dimension^2)
y_vals = reform(y_vals, obs.dimension^2)
;transform to ra/dec, back to orthoslant. For different observations obs would be different
;between the two
apply_astrometry, obs, x_arr=x_vals, y_arr=y_vals, ra_arr=ra_vals, dec_arr=dec_vals, /xy2ad
apply_astrometry, obs, ra_arr=ra_vals, dec_arr=dec_vals, x_arr=xv_hpx, y_arr=yv_hpx, /ad2xy
Comment on lines -63 to +73
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't believe you want to change the code inside this IF statement. Restore the old code, then add your new code below.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For manipulations of the obs structure you should use fhd_struct_update_obs.pro in FHD/fhd_core/setup_metadata

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To generate consistent arrays of x and y pixel values, please use:
x_vals = meshgrid(obs.dimension, obs.elements, 1)
y_vals = meshgrid(obs.dimension, obs.elements, 2)



ENDIF ELSE BEGIN
ang2vec,obs.obsdec,obs.obsra,cen_coords,/astro
Query_disc,nside,cen_coords,radius,hpx_inds0,ninds,/deg
Expand Down Expand Up @@ -90,27 +100,29 @@ ENDIF ELSE BEGIN
ENDELSE

; Test for pixels past the horizon. We don't need to be precise with this, so turn off precession, etc..
Eq2Hor,pix_ra, pix_dec, obs.JD0, alt_arr, az_arr, nutate=0,precess=0,aberration=0, refract=0, lon=obs.lon, alt=obs.alt, lat=obs.lat
horizon_i = where(alt_arr LE 0, n_horizon, complement=h_use)
IF n_horizon GT 0 THEN BEGIN
print,String(format='("Cutting ",A, " HEALPix pixels that were below the horizon.")',Strn(n_horizon))
xv_hpx = xv_hpx[h_use]
yv_hpx = yv_hpx[h_use]
hpx_inds = hpx_inds[h_use]
ENDIF
;Eq2Hor,pix_ra, pix_dec, obs.JD0, alt_arr, az_arr, nutate=0,precess=0,aberration=0, refract=0, lon=obs.lon, alt=obs.alt, lat=obs.lat
;horizon_i = where(alt_arr LE 0, n_horizon, complement=h_use)
;IF n_horizon GT 0 THEN BEGIN
; print, "im here"
; print,String(format='("Cutting ",A, " HEALPix pixels that were below the horizon.")',Strn(n_horizon))
; xv_hpx = xv_hpx[h_use]
; yv_hpx = yv_hpx[h_use]
; hpx_inds = hpx_inds[h_use]
;ENDIF

x_frac=1.-(xv_hpx-Floor(xv_hpx))
y_frac=1.-(yv_hpx-Floor(yv_hpx))

min_bin=Min(Floor(xv_hpx)+dimension*Floor(yv_hpx))>0L
max_bin=Max(Ceil(xv_hpx)+dimension*Ceil(yv_hpx))<(dimension*elements-1L)
max_bin=Max(Ceil(xv_hpx)+dimension*Ceil(yv_hpx))<(dimension*elements-1L)
; give hist some data, min and max for bins, binsize=1?
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The IDL histogram is a very powerful tool, but is also one of the most confusing. I strongly recommend reading http://www.idlcoyote.com/tips/histogram_tutorial.html

h00=histogram(Floor(xv_hpx)+dimension*Floor(yv_hpx),min=min_bin,max=max_bin,/binsize,reverse_ind=ri00)
h01=histogram(Floor(xv_hpx)+dimension*Ceil(yv_hpx),min=min_bin,max=max_bin,/binsize,reverse_ind=ri01)
h10=histogram(Ceil(xv_hpx)+dimension*Floor(yv_hpx),min=min_bin,max=max_bin,/binsize,reverse_ind=ri10)
h11=histogram(Ceil(xv_hpx)+dimension*Ceil(yv_hpx),min=min_bin,max=max_bin,/binsize,reverse_ind=ri11)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The above code is calculating the four corners surrounding every pixel, and putting the results in difficult-to-understand lookup tables (the "reverse indices").

htot=h00+h01+h10+h11
inds=where(htot,n_img_use)

n_arr=htot[inds]

i_use=inds+min_bin
Expand Down Expand Up @@ -151,7 +163,9 @@ FOR i=0L,n_img_use-1L DO BEGIN
*ija[i]=ija0

ENDFOR

;sa=weighting function?
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ija and sa are written in the sparse matrix notation of sprsax2.pro, in FH/fhd_core/gridding . The ijas are the coordinates with data in the sparse matrix, and the sas are the values.

;set hpx inds to be the right thing for the orthoslant
hpx_inds = findgen(obs.dimension^2)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do not use findgen() to generate array indices, ever. For larger arrays it lacks sufficient precision and you will get strange results. Use Lindgen() instead.

hpx_cnv={nside:nside,ija:ija,sa:sa,i_use:i_use,inds:hpx_inds}
IF tag_exist(obs,'healpix') THEN BEGIN
IF N_Elements(restrict_hpx_inds) NE 1 THEN ind_list="UNSPECIFIED" ELSE ind_list=restrict_hpx_inds
Expand Down
56 changes: 49 additions & 7 deletions fhd_core/HEALPix/healpix_snapshot_cube_generate.pro
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@ PRO healpix_snapshot_cube_generate,obs_in,status_str,psf_in,cal,params,vis_arr,v
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,$
restrict_hpx_inds=restrict_hpx_inds,hpx_radius=hpx_radius,cmd_args=cmd_args,save_uvf=save_uvf,save_imagecube=save_imagecube,$
obs_out=obs_out,psf_out=psf_out,ps_tile_flag_list=ps_tile_flag_list,_Extra=extra
obs_out=obs_out,psf_out=psf_out,ps_tile_flag_list=ps_tile_flag_list,beam_ptr=beam_ptr,_Extra=extra

t0=Systime(1)

Expand Down Expand Up @@ -184,31 +184,61 @@ PRO healpix_snapshot_cube_generate,obs_in,status_str,psf_in,cal,params,vis_arr,v
undefine_fhd,weights_arr1,variance_arr1,residual_arr1,dirty_arr1,model_arr1 ;free memory for beam_arr later!
IF iter EQ n_iter-1 THEN undefine_fhd,beam_arr

;set up the pointers for later arrays

dirty_img_arr=Ptrarr(n_pol,n_freq,/allocate)
weights_img_arr=Ptrarr(n_pol,n_freq,/allocate)
variance_img_arr=Ptrarr(n_pol,n_freq,/allocate)
model_img_arr=Ptrarr(n_pol,n_freq,/allocate)

dirty_uv_arr=Ptrarr(n_pol,n_freq,/allocate)
weights_uv_arr=Ptrarr(n_pol,n_freq,/allocate)
variance_uv_arr=Ptrarr(n_pol,n_freq,/allocate)
model_uv_arr=Ptrarr(n_pol,n_freq,/allocate)

FOR pol_i=0,n_pol-1 DO BEGIN
IF dirty_flag THEN BEGIN
dirty_cube=fltarr(n_hpx,n_freq_use)
;write index in much more efficient memory access order
FOR fi=Long64(0),n_freq_use-1 DO dirty_cube[n_hpx*fi]=Temporary(*dirty_hpx_arr[pol_i,fi])
;remove temporary to avoid losing arrays
FOR fi=Long64(0),n_freq_use-1 DO dirty_cube[n_hpx*fi]=*dirty_hpx_arr[pol_i,fi]
ENDIF

IF model_flag THEN BEGIN
model_cube=fltarr(n_hpx,n_freq_use)
FOR fi=Long64(0),n_freq_use-1 DO model_cube[n_hpx*fi]=Temporary(*model_hpx_arr[pol_i,fi])
FOR fi=Long64(0),n_freq_use-1 DO model_cube[n_hpx*fi]=*model_hpx_arr[pol_i,fi]
ENDIF

IF residual_flag THEN BEGIN
res_cube=fltarr(n_hpx,n_freq_use)
FOR fi=Long64(0),n_freq_use-1 DO res_cube[n_hpx*fi]=Temporary(*residual_hpx_arr[pol_i,fi])
FOR fi=Long64(0),n_freq_use-1 DO res_cube[n_hpx*fi]=*residual_hpx_arr[pol_i,fi]
ENDIF

weights_cube=fltarr(n_hpx,n_freq_use)
FOR fi=Long64(0),n_freq_use-1 DO weights_cube[n_hpx*fi]=Temporary(*weights_hpx_arr[pol_i,fi])
FOR fi=Long64(0),n_freq_use-1 DO weights_cube[n_hpx*fi]=*weights_hpx_arr[pol_i,fi]

variance_cube=fltarr(n_hpx,n_freq_use)
FOR fi=Long64(0),n_freq_use-1 DO variance_cube[n_hpx*fi]=Temporary(*variance_hpx_arr[pol_i,fi])
FOR fi=Long64(0),n_freq_use-1 DO variance_cube[n_hpx*fi]=*variance_hpx_arr[pol_i,fi]

beam_squared_cube=fltarr(n_hpx,n_freq_use)
FOR fi=Long64(0),n_freq_use-1 DO beam_squared_cube[n_hpx*fi]=Temporary(*beam_hpx_arr[pol_i,fi])
FOR fi=Long64(0),n_freq_use-1 DO beam_squared_cube[n_hpx*fi]=*beam_hpx_arr[pol_i,fi]

;now do the FFTs back and normalize

FOR fi=Long64(0),n_freq_use-1 DO *dirty_img_arr[pol_i,fi]=reform(*dirty_hpx_arr[pol_i,fi], SQRT(n_hpx), SQRT(n_hpx))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't use SQRT(n_hpx). Instead use obs.dimension, obs.elements

FOR fi=Long64(0),n_freq_use-1 DO *dirty_uv_arr[pol_i,fi]=fft_shift(fft(fft_shift((*dirty_img_arr[pol_i,fi])*(obs_out.degpix*!DtoR)^2), /inverse))


FOR fi=Long64(0),n_freq_use-1 DO *weights_img_arr[pol_i,fi]=reform(*weights_hpx_arr[pol_i,fi], SQRT(n_hpx), SQRT(n_hpx))
FOR fi=Long64(0),n_freq_use-1 DO *weights_uv_arr[pol_i,fi]=fft_shift(fft(fft_shift((*weights_img_arr[pol_i,fi])*(obs_out.degpix*!DtoR)^2), /inverse))


FOR fi=Long64(0),n_freq_use-1 DO *variance_img_arr[pol_i,fi]=reform(*variance_hpx_arr[pol_i,fi], SQRT(n_hpx), SQRT(n_hpx))
FOR fi=Long64(0),n_freq_use-1 DO *variance_uv_arr[pol_i,fi]=fft_shift(fft(fft_shift((*variance_img_arr[pol_i,fi])*(obs_out.degpix*!DtoR)^2), /inverse))

FOR fi=Long64(0),n_freq_use-1 DO *model_img_arr[pol_i,fi]=reform(*model_hpx_arr[pol_i,fi], SQRT(n_hpx), SQRT(n_hpx))
FOR fi=Long64(0),n_freq_use-1 DO *model_uv_arr[pol_i,fi]=fft_shift(fft(fft_shift((*model_img_arr[pol_i,fi])*(obs_out.degpix*!DtoR)^2), /inverse))


;call fhd_save_io first to obtain the correct path. Will NOT update status structure yet
fhd_save_io,status_str,file_path_fhd=file_path_fhd,var=cube_name[iter],pol_i=pol_i,path_use=path_use,/no_save,_Extra=extra
Expand All @@ -219,9 +249,21 @@ 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
dirty_cube=(model_cube=(res_cube=(weights_cube=(variance_cube=(beam_squared_cube=0)))))
ENDFOR

image_filepath = file_path_fhd+ '_' + uvf_name[iter] +'_interped_images.sav'
save, filename = image_filepath, dirty_img_arr, weights_img_arr, variance_img_arr, model_img_arr, obs_out, /compress

uvf_filepath = file_path_fhd+ '_' + uvf_name[iter] +'_gridded_uvf.sav'
save, filename = uvf_filepath, dirty_uv_arr, weights_uv_arr, variance_uv_arr, model_uv_arr, obs_out, /compress

undefine_fhd,dirty_hpx_arr,model_hpx_arr,weights_hpx_arr,variance_hpx_arr,beam_hpx_arr,residual_hpx_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
Expand Down
9 changes: 9 additions & 0 deletions fhd_utils/FFT/fft_shift_cube.pro
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
FUNCTION fft_shift_cube, image_cube
dimension=(size(image_cube,/dimension))[0]
elements=(size(image_cube,/dimension))[1]
n_slice=(size(image_cube,/dimension))[2]
shifted_image_cube = Fltarr(dimension, elements, n_slice)
FOR s=0,n_slice-1 DO $
shifted_image_cube[*, *, s] = fft_shift(Reform(image_cube[*, *, s]))
RETURN, shifted_image_cube
END
Empty file added idl_startup.pro
View file
Open in desktop
Empty file.