Reducing skymodel confusion

dictionary.md

@@ -225,7 +225,6 @@ This is a work in progress; please add keywords as you find them in alphabetical
   -*Default*: equal to `max_baseline` <br />
 
 **max_calibration_sources**: limits the number of sources used in the calibration. Sources are weighted by apparent brightness before applying the cut. Note that extended sources with many associated source components count as only a single source. <br />
-  -*Dependency*: The sources are also included in the model if `return_cal_visibilities` is set. <br />
   -*Default*: All valid sources in the catalog are used. <br />
 
 **min_cal_baseline**: the minimum baseline length in wavelengths to be used in calibration. <br />
@@ -417,6 +416,11 @@ WARNING! Options in this section may change without notice, and should never be
   -*Default*: Unset, unless model_visibilities and model_subtract_sidelobe_catalog are set OR deconvolve and subtract_sidelobe_catalog are set <br />
   -*eor_wrapper_defaults*: 1 <br />
 
+**combine_skymodels**: combines the calibration and model skymodel together. Sources are marked as duplicates if matched by 1/100th of a pixel. If `return_cal_visibilities` is set, then only the non-duplicate model sources are calculated, which are added to the returned calibration visibilities to build the full model. This assumes that the calibration skymodel is a subset of the model skymodel. This reduces the number of calculations and is thus more efficient. If `return_cal_visibilities` is not set, then the entire model skymodel is calculated. If `save_skymodel` is set, then the calibration skymodel (marked with `include_calibration=1` in the structure) and non-duplicate model skymodel (marked with `include_calibration=0` in the structure) are saved as one output. <br />
+  -*Dependency*: `calibrate_visibilities` and `model_visibilities` must be set to 1 <br />
+  -*Default*: Unset <br />
+  -*Turn off/on*: 0/1 <br />
+
 **conserve_memory**: split the model DFT into matrix chunks to reduce the memory load at the cost of extra walltime. If set to greater than 1e6, then it will set the maximum memory used in bytes. <br />
   -*Default*: 1 (100 Mb) <br />
   -*Turn off/on*: 0/1 (optionally set to bytes)<br />
@@ -441,7 +445,7 @@ WARNING! Options in this section may change without notice, and should never be
   -*Default*: 0 <br />
 
 **max_model_sources**: limits the number of sources used in the model. Sources are weighted by apparent brightness before applying the cut. Note that extended sources with many associated source components count as only a single source. <br />
-  -*Dependency*: `model_visibilities` must be set to 1 in order for the keyword to take effect. If `return_cal_visibilities` is set, then the final model will include all calibration sources and all model sources (duplicates are caught and included only once). <br />
+  -*Dependency*: `model_visibilities` must be set to 1 in order for the keyword to take effect. <br />
   -*Default*: All valid sources are used. !Q <br />
 
 **model_catalog_file_path**: a catalog of sources to be used to make model visibilities for subtraction. The source catalog can be provided as either a IDL .sav file or a .skyh5 file. A .skyh5 file must be created with pyradiosky or conform to the pyradiosky formatting convention.<br />
@@ -510,6 +514,10 @@ WARNING! Options in this section may change without notice, and should never be
   -*Default*: 0 <br />
   -*eor_wrapper_defaults*: 10.x`pad_uv_image` <br />
 
+**save_skymodel**: save the skymodel. If both `calibrate_visibilities` and `model_visibilities` are set, then save them individually as `skymodel_cal` and `skymodel_model` unless `combine_skymodels` is set. <br />
+  -*Turn off/on*: 0/1 <br />
+  -*Default*: 0 <br />
+
 **save_uvf**: saves the gridded uv plane as a function of frequency for dirty, model, weights, and variance cubes. <br />
   -*Turn off/on*: 0/1 <br />
   -*Default*: 0 <br />

fhd_core/fhd_main.pro

@@ -5,7 +5,7 @@ PRO fhd_main, file_path_vis, status_str, export_images=export_images, cleanup=cl
     file_path_fhd=file_path_fhd, force_data=force_data, force_no_data=force_no_data, freq_start=freq_start, freq_end=freq_end,$
     calibrate_visibilities=calibrate_visibilities, transfer_calibration=transfer_calibration, error=error,$
     calibration_catalog_file_path=calibration_catalog_file_path, dft_threshold=dft_threshold,$
-    calibration_visibilities_subtract=calibration_visibilities_subtract,$
+    calibration_visibilities_subtract=calibration_visibilities_subtract,save_skymodel=save_skymodel,combine_skymodels=combine_skymodels,$
     weights_grid=weights_grid, save_visibilities=save_visibilities, return_cal_visibilities=return_cal_visibilities,$
     return_decon_visibilities=return_decon_visibilities, snapshot_healpix_export=snapshot_healpix_export, cmd_args=cmd_args, log_store=log_store,$
     generate_vis_savefile=generate_vis_savefile, model_visibilities=model_visibilities, model_catalog_file_path=model_catalog_file_path,$
@@ -100,7 +100,7 @@ IF data_flag LE 0 THEN BEGIN
         IF ~Keyword_Set(calibration_source_list) THEN $
           calibration_source_list=generate_source_cal_list(obs,psf,catalog_path=catalog_use,_Extra=extra)
       ENDIF
-      skymodel_cal=fhd_struct_init_skymodel(obs,source_list=calibration_source_list,catalog_path=catalog_use,return_cal=1,diffuse_model=diffuse_calibrate,_Extra=extra)
+      skymodel_cal=fhd_struct_init_skymodel(obs,source_list=calibration_source_list,catalog_path=catalog_use,calibration_flag=1,diffuse_model=diffuse_calibrate,_Extra=extra)
       cal=fhd_struct_init_cal(obs,params,skymodel_cal,source_list=calibration_source_list,$
         catalog_path=catalog_use,transfer_calibration=transfer_calibration,_Extra=extra)
     ENDIF
@@ -175,20 +175,49 @@ IF data_flag LE 0 THEN BEGIN
         fhd_save_io,status_str,vis_weights,var='vis_weights',/compress,file_path_fhd=file_path_fhd,_Extra=extra
 
     IF Keyword_Set(model_visibilities) THEN BEGIN
-        IF Keyword_Set(model_catalog_file_path) AND ~Keyword_set(model_uv_transfer) THEN BEGIN
-            model_source_list=generate_source_cal_list(obs,psf,catalog_path=model_catalog_file_path,/model_visibilities,_Extra=extra) 
-            IF Keyword_Set(return_cal_visibilities) OR Keyword_Set(calibration_visibilities_subtract) THEN $
+        print, "Modelling visibilities separate from visibilities used for calibration"
+        IF ~Keyword_set(model_uv_transfer) THEN BEGIN
+            ;Build the source list for modelling visilibilities
+            model_source_list=generate_source_cal_list(obs,psf,jones=jones,catalog_path=model_catalog_file_path,/model_visibilities,_Extra=extra) 
+            IF Keyword_Set(return_cal_visibilities) AND Keyword_Set(combine_skymodels) THEN BEGIN
+                ;Combine the calibration source list with the model source list if calibration visibilities still exist
+                ;and only calculate excess model source visibilities and add them to the calibration visibilities.
+                ;Great for avoiding excess calculations with large, semi-matched input catalogs.
                 model_source_list=source_list_append(obs,model_source_list,skymodel_cal.source_list,/exclude)
+                skymodel_model=fhd_struct_init_skymodel(obs,source_list=model_source_list,catalog_path=model_catalog_file_path,$
+                    diffuse_model=diffuse_model,return_cal=return_cal_visibilities,_Extra=extra)
+                vis_model_arr=vis_source_model(skymodel_model,obs,status_str,psf,params,vis_weights,0,jones,model_uv_arr=model_uv_arr,$
+                    timing=model_timing,silent=silent,error=error,vis_model_ptr=vis_model_arr,calibration_flag=0,file_path_fhd=file_path_fhd,_Extra=extra)
+            ENDIF ELSE BEGIN
+                ;Build the skymodel and model visilibilities for all input model sources.
+                skymodel_model=fhd_struct_init_skymodel(obs,source_list=model_source_list,catalog_path=model_catalog_file_path,$
+                    diffuse_model=diffuse_model,return_cal=return_cal_visibilities,_Extra=extra)
+                vis_model_arr=vis_source_model(skymodel_model,obs,status_str,psf,params,vis_weights,0,jones,model_uv_arr=model_uv_arr,$
+                    timing=model_timing,silent=silent,error=error,calibration_flag=0,file_path_fhd=file_path_fhd,_Extra=extra)
+            ENDELSE
         ENDIF
-        skymodel_update=fhd_struct_init_skymodel(obs,source_list=model_source_list,catalog_path=model_catalog_file_path,$
-            diffuse_model=diffuse_model,return_cal=return_cal_visibilities,_Extra=extra)
-        vis_model_arr=vis_source_model(skymodel_update,obs,status_str,psf,params,vis_weights,0,jones,model_uv_arr=model_uv_arr,$
-            timing=model_timing,silent=silent,error=error,vis_model_ptr=vis_model_arr,calibration_flag=0,file_path_fhd=file_path_fhd,_Extra=extra)
     ENDIF
 
-    IF Keyword_Set(skymodel_cal) OR Keyword_Set(skymodel_update) THEN $
-      skymodel=fhd_struct_combine_skymodel(obs, skymodel_cal, skymodel_update,calibration_flag=(Keyword_Set(return_cal_visibilities) OR Keyword_Set(calibration_visibilities_subtract)))
-    fhd_save_io,status_str,skymodel,var='skymodel',/compress,file_path_fhd=file_path_fhd,_Extra=extra
+    IF Keyword_Set(skymodel_cal) OR Keyword_Set(skymodel_model) THEN BEGIN
+        ;Combine and/or save the skymodels if desired.
+        ;If there is one skymodel, it will be saved as 'skymodel'. Otherwise, they will be saved as 'skymodel_cal' and 'skymodel_model'.
+        IF Keyword_Set(combine_skymodels) THEN BEGIN
+            skymodel=fhd_struct_combine_skymodel(obs, skymodel_cal, skymodel_model)
+            IF Keyword_set(save_skymodel) THEN fhd_save_io,status_str,skymodel,var='skymodel',/compress,file_path_fhd=file_path_fhd,_Extra=extra
+        ENDIF ELSE BEGIN
+            IF Keyword_set(save_skymodel) THEN BEGIN
+                IF Keyword_Set(skymodel_cal) AND ~Keyword_Set(skymodel_model) THEN BEGIN
+                    skymodel = temporary(skymodel_cal)
+                    fhd_save_io,status_str,skymodel,var='skymodel',/compress,file_path_fhd=file_path_fhd,_Extra=extra
+                ENDIF ELSE BEGIN
+                    fhd_save_io,status_str,skymodel_cal,var='skymodel_cal',/compress,file_path_fhd=file_path_fhd,_Extra=extra
+                    fhd_save_io,status_str,skymodel_model,var='skymodel_model',/compress,file_path_fhd=file_path_fhd,_Extra=extra
+                    ;Rename one of the skymodels to pass along keywords to deconvolution
+                    skymodel = temporary(skymodel_cal)
+                ENDELSE
+            ENDIF
+        ENDELSE
+    ENDIF
 
     IF N_Elements(vis_model_arr) LT n_pol THEN vis_model_arr=Ptrarr(n_pol) ;supply as array of null pointers to allow it to be indexed, but not be used
     model_flag=min(Ptr_valid(vis_model_arr))
@@ -281,7 +310,7 @@ IF Keyword_Set(production) THEN BEGIN
 ENDIF
 
 undefine_fhd,map_fn_arr,image_uv_arr,weights_arr,model_uv_arr,vis_arr,vis_weights,vis_model_arr
-undefine_fhd,obs,cal,jones,layout,psf,antenna,fhd_params,skymodel,skymodel_cal,skymodel_update
+undefine_fhd,obs,cal,jones,layout,psf,antenna,fhd_params,skymodel,skymodel_cal,skymodel_model
 
 run_report, start_mem, t0, silent=silent
 print,''

fhd_core/setup_metadata/fhd_save_io.pro

@@ -50,7 +50,7 @@ IF Keyword_Set(text) THEN BEGIN
 ENDIF
 
 IF N_Elements(var_name) EQ 0 THEN var_name='' ELSE var_name=StrLowCase(var_name)
-IF var_name NE 'status_str' THEN IF not tag_exist(status_use,var_name) THEN RETURN
+IF var_name NE 'status_str' THEN IF ~strmatch(var_name,'skymodel*') THEN IF not tag_exist(status_use,var_name) THEN RETURN
 
 CASE var_name OF ;listed in order typically generated
     'status_str':BEGIN path_add='_status' & subdir='metadata' & END
@@ -63,6 +63,8 @@ CASE var_name OF ;listed in order typically generated
     'jones':BEGIN status_use.jones=1 & path_add='_jones' & subdir='beams'& END
     'cal':BEGIN status_use.cal=1 & path_add='_cal' & subdir='calibration'& END
     'skymodel':BEGIN status_use.skymodel=1 & path_add='_skymodel' & subdir='output_data'& END
+    'skymodel_cal':BEGIN status_use.skymodel=1 & path_add='_skymodel_cal' & subdir='output_data'& END
+    'skymodel_model':BEGIN status_use.skymodel=1 & path_add='_skymodel_model' & subdir='output_data'& END
     'source_array':BEGIN status_use.source_array=1 & path_add='_source_array' & subdir='output_data'& END
     'vis_weights':BEGIN
         IF Tag_exist(status_use,"vis_weights") THEN status_use.vis_weights=1 ELSE status_use.flag_arr=1

fhd_core/source_modeling/fhd_struct_combine_skymodel.pro

@@ -1,8 +1,16 @@
-FUNCTION fhd_struct_combine_skymodel, obs, skymodel_cal, skymodel_update, calibration_flag=calibration_flag
+FUNCTION fhd_struct_combine_skymodel, obs, skymodel_cal, skymodel_update
 
-IF N_Elements(skymodel_cal) EQ 0 THEN RETURN,skymodel_update
-IF N_Elements(skymodel_update) EQ 0 THEN RETURN,skymodel_cal
-IF N_Elements(calibration_flag) EQ 0 THEN calibration_flag=1 ;set to indicate skymodel_cal really is the calibration model
+IF N_Elements(skymodel_cal) EQ 0 THEN RETURN, skymodel_update
+IF N_Elements(skymodel_update) EQ 0 THEN RETURN, skymodel_cal
+
+;Combine calibration flag arrays to show which sources were used to calibrate
+if N_elements(skymodel_cal.include_calibration) EQ 1 THEN BEGIN
+    cal_flag = INTARR(skymodel_cal.n_sources) + skymodel_cal.include_calibration
+ENDIF ELSE cal_flag = skymodel_cal.include_calibration
+if N_elements(skymodel_update.include_calibration) EQ 1 THEN BEGIN
+    model_flag = INTARR(skymodel_update.n_sources) + skymodel_update.include_calibration
+ENDIF ELSE model_flag = skymodel_update.include_calibration
+calibration_flag = [cal_flag, model_flag]
 
 source_list=source_list_append(obs,skymodel_update.source_list,skymodel_cal.source_list)
 IF Keyword_Set(skymodel_update.diffuse_model) THEN BEGIN
@@ -26,7 +34,7 @@ ENDIF ELSE BEGIN
 ENDELSE
 
 skymodel=fhd_struct_init_skymodel(obs,source_list=source_list,$
-    return_cal=calibration_flag,catalog_path=catalog_path,$
+    calibration_flag=calibration_flag,catalog_path=catalog_path,$
     galaxy_model=galaxy_model,galaxy_spectral_index=galaxy_spectral_index,$
     diffuse_model=diffuse_model,diffuse_spectral_index=diffuse_spectral_index)
 

fhd_core/source_modeling/fhd_struct_init_skymodel.pro

@@ -1,9 +1,9 @@
 FUNCTION fhd_struct_init_skymodel,obs,source_list=source_list,catalog_path=catalog_path,$
     galaxy_model=galaxy_model,galaxy_spectral_index=galaxy_spectral_index,$
     diffuse_model=diffuse_model,diffuse_spectral_index=diffuse_spectral_index,$
-    return_cal_visibilities=return_cal_visibilities,_Extra=extra
+    calibration_flag=calibration_flag,_Extra=extra
 
-IF N_Elements(return_cal_visibilities) EQ 0 THEN calibration_flag=0 ELSE calibration_flag=return_cal_visibilities
+IF N_Elements(calibration_flag) EQ 0 THEN calibration_flag=0
 IF Keyword_Set(galaxy_model) THEN galaxy_flag=1 ELSE galaxy_flag=0
 IF Keyword_Set(diffuse_model) THEN diffuse_filepath=diffuse_model ELSE diffuse_filepath=''
 IF Keyword_Set(catalog_path) THEN catalog_path_use=file_basename(catalog_path,'.sav',/fold_case) ELSE catalog_path_use=''

fhd_core/source_modeling/source_list_append.pro

@@ -15,14 +15,11 @@ ENDFOR
 
 dup_i=where(dup_flag GE 0,n_dup,complement=uniq_i1,ncomplement=n_uniq1)
 IF Keyword_Set(exclude_duplicate) THEN BEGIN
-    IF n_uniq1 GT 0 THEN source_list=source_list1[uniq_i1] ELSE source_list=source_comp_init(n_sources=0)
+    ;Only use unique sources from first source list.
+    ;If no unique sources, return the entire first source list.
+    IF n_uniq1 GT 0 THEN source_list=source_list1[uniq_i1] ELSE source_list=source_list1
 ENDIF ELSE BEGIN
-    CASE n_dup OF
-        0:source_list=[source_list1,source_list2] 
-        n2:source_list=source_list1
-        n1:source_list=source_list2
-        ELSE:source_list=[source_list1[uniq_i1],source_list2]
-    ENDCASE
+    source_list=[source_list1,source_list2] 
 ENDELSE
 
 RETURN,source_list