Merge pull request #13 from miguelcarcamov/development

Development

pyproject.toml

@@ -13,7 +13,7 @@ maintainers = [
 description = "A Python object oriented framework to do radio-interferometric self-calibration"
 readme = "README.md"
 license = {text = "GNU GPL"}
-requires-python = ">=3.8, <=3.8.10"
+requires-python = ">=3.8, <=3.8.13"
 classifiers = [
     "Programming Language :: Python :: 3",
     "Operating System :: OS Independent"

requirements.txt

@@ -1,5 +1,5 @@
 almatasks==1.5.2
-astropy==5.1
+astropy==5.1.1
 casadata==2022.9.5
 casampi==0.5.01
 casaplotms==1.8.7
@@ -10,8 +10,8 @@ casatestutils==6.5.2.26
 casatools==6.5.2.26
 casaviewer==1.6.6
 mpi4py==3.1.3
-numpy==1.23
-git+https://gitlab.com/miguelcarcamov/pyralysis@development#egg=pyralysis
+numpy==1.23.4
+git+https://gitlab.com/clirai/pyralysis@development#egg=pyralysis
 python-casacore==3.5.2
 reproject==0.9.1
 scipy==1.9.1

src/snow/__init__.py

@@ -0,0 +1,10 @@
+try:
+    # -- Distribution mode --
+    # import from _version.py generated by setuptools_scm during release
+    from ._version import version as __version__
+except ImportError:
+    # -- Source mode --
+    # use setuptools_scm to get the current version from src using git
+    from setuptools_scm import get_version as _gv
+    from os import path as _path
+    __version__ = _gv(_path.join(_path.dirname(__file__), _path.pardir))

src/snow/imaging/__init__.py

@@ -1,4 +1,4 @@
-from .gpuvmem import *
-from .imager import *
-from .tclean import *
-from .wsclean import *
+from .gpuvmem import GPUvmem
+from .imager import Imager
+from .tclean import Tclean
+from .wsclean import WSClean

src/snow/imaging/gpuvmem.py

@@ -328,12 +328,15 @@ def __restore(self,
 
         importfits(imagename="model_out", fitsimage=model_fits, overwrite=True)
 
+        aux_reference_freq = self._check_reference_frequency()
+
         tclean(
             vis=residual_ms,
             imagename=residual_image,
             specmode='mfs',
             deconvolver='hogbom',
             niter=0,
+            reffreq=aux_reference_freq,
             stokes=self.stokes,
             nterms=1,
             weighting=self.weighting,
@@ -409,12 +412,14 @@ def __create_model_input(self, name="model_input") -> str:
             A string with absolute path to the output FITS image file
         """
         fits_image = name + '.fits'
+        aux_reference_freq = self._check_reference_frequency()
         tclean(
             vis=self.inputvis,
             imagename=name,
             specmode='mfs',
             niter=0,
             deconvolver='hogbom',
+            reffreq=aux_reference_freq,
             interactive=False,
             cell=self.cell,
             stokes=self.stokes,
@@ -436,7 +441,15 @@ def run(self, imagename=""):
             Absolute path to the output image name file
         """
         if self.model_input is None:
+            print("Model input is None - Creating model_input...")
             self.model_input = self.__create_model_input(imagename + "_input")
+
+        if self.reference_freq is None:
+            print("Reference frequency is None - Reading CRVAL3 as reference frequency...")
+            with fits.open(self.model_input) as hdul:
+                header = hdul[0].header
+                self.reference_freq = Quantity(header['CRVAL3'] * u.Hz)
+
         model_output = imagename + ".fits"
         _residual_output = imagename + "_" + self.residual_output
         restored_image = imagename + ".restored"
@@ -457,6 +470,18 @@ def run(self, imagename=""):
         if self.gridding:
             args += " -g " + str(self.gridding_threads)
 
+        if self.reference_freq:
+            if isinstance(self.reference_freq, Quantity):
+                args += " -F " + str(self.reference_freq.to(u.Hz).value)
+            elif isinstance(self.reference_freq, float):
+                args += " -F " + str(self.reference_freq)
+            elif isinstance(self.reference_freq, str):
+                args += " -F " + self.reference_freq
+            else:
+                raise NotImplementedError(
+                    "Type {} has not implementation in snow".format(type(self.reference_freq))
+                )
+
         if self.print_images:
             args += " --print-images"
 

src/snow/imaging/imager.py

@@ -2,6 +2,10 @@
 from dataclasses import dataclass
 from dataclasses import field as _field
 
+from typing import Union
+
+from astropy.units import Quantity
+
 from ..utils import (calculate_number_antennas, calculate_psnr_fits, calculate_psnr_ms)
 
 
@@ -62,6 +66,7 @@ class Imager(metaclass=ABCMeta):
     data_column: str = "corrected"
     M: int = 512
     N: int = 512
+    reference_freq: Union[str, float, Quantity, None] = None
     niter: int = 100
     noise_pixels: int = None
     save_model: bool = True
@@ -128,6 +133,21 @@ def _calculate_statistics_msimage(
         self.peak = peak
         self.stdv = stdv
 
+    def _check_reference_frequency(self):
+        aux_reference_freq = ""
+        if self.reference_freq:
+            if isinstance(self.reference_freq, Quantity) or isinstance(self.reference_freq, float):
+                aux_reference_freq = str(self.reference_freq)
+            elif self.reference_freq is None:
+                aux_reference_freq = ""
+            elif isinstance(self.reference_freq, str):
+                aux_reference_freq = self.reference_freq
+            else:
+                raise NotImplementedError(
+                    "Type {} has not implementation in snow".format(type(self.reference_freq))
+                )
+        return aux_reference_freq
+
     @abstractmethod
     def run(self, imagename=""):
         return

src/snow/imaging/tclean.py

@@ -75,6 +75,7 @@ class Tclean(Imager):
     noise_threshold: float = 4.25
     sidelobe_threshold: float = 2.0
     min_beam_frac: float = 0.3
+    grow_iterations: float = 75.0
     specmode: str = ""
     gridder: str = "standard"
     wproj_planes: int = -1
@@ -95,6 +96,7 @@ def __post_init__(self):
 
     def run(self, imagename=""):
         __imsize = [self.M, self.N]
+        aux_reference_freq = self._check_reference_frequency()
         tclean(
             vis=self.inputvis,
             imagename=imagename,
@@ -108,6 +110,7 @@ def run(self, imagename=""):
             scales=self.scales,
             nterms=self.nterms,
             imsize=__imsize,
+            reffreq=aux_reference_freq,
             cell=self.cell,
             weighting=self.weighting,
             robust=self.robust,
@@ -120,12 +123,13 @@ def run(self, imagename=""):
             pbcor=self.pbcor,
             uvtaper=self.uvtaper,
             savemodel=self.clean_savemodel,
-            usemask=self.usemask,
+            usemask=self.use_mask,
             negativethreshold=self.negative_threshold,
             lownoisethreshold=self.low_noise_threshold,
             noisethreshold=self.noise_threshold,
             sidelobethreshold=self.sidelobe_threshold,
             minbeamfrac=self.min_beam_frac,
+            growiterations=self.grow_iterations,
             cycleniter=self.cycle_niter,
             verbose=self.verbose
         )

src/snow/selfcalibration/__init__.py

@@ -1,4 +1,4 @@
-from .ampcal import *
-from .apcal import *
-from .phasecal import *
-from .selfcal import *
+from .ampcal import Ampcal
+from .apcal import AmpPhasecal
+from .phasecal import Phasecal
+from .selfcal import Selfcal

src/snow/selfcalibration/selfcal.py

@@ -388,7 +388,7 @@ def _finish_selfcal_iteration(self, current_iteration: int = 0) -> bool:
             return False
 
     def _flag_dataset(
-        self, datacolumn="RESIDUAL", mode="rflag", timedevscale=3.0, freqdevscale=3.0
+        self, datacolumn=None, mode="rflag", timedevscale=3.0, freqdevscale=3.0
     ) -> None:
         """
         Protected method that flag the dataset on each iteration. This functions aims to flag residual outliers.
@@ -411,16 +411,23 @@ def _flag_dataset(
             For spectral analysis, flag a point if local rms around it is larger than freqdevscale $x$ freqdev.
         """
 
+        if datacolumn is None:
+            if is_column_in_ms(self.visfile, "CORRECTED_DATA"):
+                datacolumn = "residual"
+            else:
+                datacolumn = "residual_data"
+
+        print("Flagging {0} data column using {1}".format(datacolumn, mode))
+
         flagdata(
             vis=self.visfile,
             mode=mode,
             datacolumn=datacolumn,
             field=self.imager.field,
             timecutoff=5.0,
             freqcutoff=5.0,
-            timefit='line',
             freqfit='line',
-            flagdimension='freqtime',
+            flagdimension='freq',
             extendflags=False,
             timedevscale=timedevscale,
             freqdevscale=freqdevscale,
@@ -429,6 +436,7 @@ def _flag_dataset(
             growaround=False,
             flagneartime=False,
             flagnearfreq=False,
+            ntime="scan",
             action='apply',
             flagbackup=True,
             overwrite=True,
@@ -494,7 +502,7 @@ def _plot_selfcal(
             # gridcols=subplot[1], antenna=antenna, timerange=timerange, plotrange=plotrange, plotfile=figfile_name,
             # overwrite=True, showgui=False)
 
-    def selfcal_output(self, overwrite=False, _statwt=False) -> str:
+    def selfcal_output(self, overwrite=False, _statwt=False, min_samp=8) -> str:
         """
         Public function that creates a new measurement set only taking the corrected column.
         If _statwt is True then applies the statwt function and creates a .statwt measurement
@@ -505,7 +513,8 @@ def selfcal_output(self, overwrite=False, _statwt=False) -> str:
             Whether to overwrite the measurement set files
         _statwt :
             Whether to create a new measurement applying the statwt function
-
+        min_samp :
+            Minimum number of unflagged visibilities for estimating the scatter if _statwt is True
         Returns
         -------
         Name of the self-calibrated measurement set file
@@ -529,7 +538,7 @@ def selfcal_output(self, overwrite=False, _statwt=False) -> str:
             if os.path.exists(statwt_path):
                 shutil.rmtree(statwt_path)
             shutil.copytree(output_vis, statwt_path)
-            statwt(vis=statwt_path, datacolumn="data")
+            statwt(vis=statwt_path, datacolumn="data", minsamp=min_samp)
         return output_vis
 
     def _uvsubtract(self):

src/snow/utils/__init__.py

@@ -1,2 +1,2 @@
-from .image_utils import *
-from .selfcal_utils import *
+from .image_utils import nanrms, rms, get_header, get_hdu, get_hdul, get_data, get_header_and_data, export_ms_to_fits, calculate_psnr_fits, calculate_psnr_ms, reproject
+from .selfcal_utils import is_column_in_ms, get_table_rows, calculate_number_antennas