Bring up to date with focalpy38

Dockerfile

@@ -1,7 +1,7 @@
-FROM ubuntu:18.04 AS spython-base
+FROM ubuntu:20.04 AS spython-base
 
 LABEL "APPLICATION_NAME_BASE"="Ubuntu LTS + MeqSilhouette + dependencies"
-LABEL "OS_VERSION"="18.04"
+LABEL "OS_VERSION"="20.04"
 LABEL "SYSTEM_NAME"="MeqSv2"
 LABEL "SYSTEM_URL"="https://github.com/rdeane/MeqSilhouette"
 LABEL "AUTHOR_NAME"="Iniyan Natarajan, Robin Hall"
@@ -15,14 +15,10 @@ ENV DEBIAN_FRONTEND=noninteractive
 
 # install utility packages and repositories
 RUN apt-get update -y
-RUN apt-get install -y wget vim python-pip gcc python unzip git time
+RUN apt-get install -y wget vim python3-pip gcc python3 unzip git time
 RUN pip uninstall -y enum34
-RUN apt-get install -y build-essential cmake gfortran g++ libncurses5-dev \
-    libreadline-dev flex bison libblas-dev liblapacke-dev libcfitsio-dev \
-    wcslib-dev libfftw3-dev libhdf5-serial-dev rsync \
-    libboost-python-dev libboost-program-options-dev libboost-program-options1.65.1 \
-    libboost-program-options1.65-dev libpython2.7-dev libxml2-dev libxslt1-dev \
-    texlive-latex-extra texlive-fonts-recommended dvipng
+RUN apt-get install -y build-essential cmake g++ rsync libboost-python-dev \
+    libboost-program-options-dev texlive-latex-extra texlive-fonts-recommended dvipng cm-super
 
 # build aatm manually
 RUN cd
@@ -31,24 +27,29 @@ RUN wget -c https://launchpad.net/aatm/trunk/0.5/+download/aatm-0.5.tar.gz
 RUN tar -xzf aatm-0.5.tar.gz
 RUN cd aatm-0.5 && ./configure && make && make install
 
-# install kern 6
+# register kern 7
 RUN apt-get install -y software-properties-common
-RUN add-apt-repository -s ppa:kernsuite/kern-6
+RUN add-apt-repository -s ppa:kernsuite/kern-7
 RUN apt-add-repository multiverse
 RUN apt-add-repository restricted
 
+# install numpy==1.21 to avoid the error "numpy has no attribute 'BitGenerator' and 'asscalar'"
+#RUN pip install numpy==1.21
+
 # install required packages
 RUN apt-get install -y \
     meqtrees \
     meqtrees-timba \
-    tigger \
     tigger-lsm \
-    python-astro-tigger \
-    python-astro-tigger-lsm \
+    python3-astro-tigger \
+    python3-astro-tigger-lsm \
     casalite \
     wsclean \
     pyxis \
-    python-casacore
+    python3-casacore
+
+# install numpy==1.21 to avoid the error "numpy has no attribute 'BitGenerator' and 'asscalar'"
+RUN pip install numpy==1.21
 
 RUN apt-get clean
 
@@ -58,11 +59,11 @@ RUN casa-config --exec update-data
 # download and install MeqSilhouette from master
 RUN cd
 RUN cd /opt
-RUN git clone --depth 1 https://github.com/rdeane/MeqSilhouette.git
+RUN git clone --branch focalpy38 --depth 1 https://github.com/rdeane/MeqSilhouette.git
 RUN cd MeqSilhouette && pip install .
 RUN cd
 
-ENV MEQTREES_CATTERY_PATH=/usr/lib/python2.7/dist-packages/Cattery
+ENV MEQTREES_CATTERY_PATH=/usr/lib/python3/dist-packages/Cattery
 ENV PATH=/usr/local/bin:$PATH
-ENV PYTHONPATH=/usr/lib/python2.7/dist-packages
+ENV PYTHONPATH=/usr/local/lib/python3.8/dist-packages:/usr/lib/python3/dist-packages
 ENV LC_ALL=C

docs/source/intro.rst

@@ -6,8 +6,8 @@ MeqSilhouette is a radio interferometry observation simulator.
 `Blecher et al. (2017) <https://ui.adsabs.harvard.edu/abs/2017MNRAS.464..143B/abstract>`_ presents an earlier version which could
 generate synthetic data and apply corruptions such as scattering by interstellar medium (ISM), atmospheric effects,
 and antenna pointing offsets. 
-The current version, `Natarajan et al. (2022) <https://ui.adsabs.harvard.edu/abs/2022MNRAS.512..490N/abstract>`_, is capable of generating fully polarised, time-variable, spectrally-resolved sky models
-and propagation path effects including (but not limited to) instrumental polarisation, updated atmospheric effects and pointing
-offsets, time-variable antenna gains, and bandpass effects.
+The current version, `Natarajan et al. (2022) <https://ui.adsabs.harvard.edu/abs/2022MNRAS.512..490N/abstract>`_, is capable of generating fully polarised,
+time-variable, spectrally-resolved sky models and propagation path effects including (but not limited to) instrumental polarisation, updated atmospheric 
+effects and pointing offsets, time-variable antenna gains, and bandpass effects.
 
 This documentation explains the steps necessary for installing and using MeqSilhouette.

docs/source/requirements.rst

@@ -2,31 +2,38 @@
 Requirements & Installation
 ===========================
 
-.. note:: This is the documentation for the default (master) branch of MeqSv2. For the Python 3 version `click here <https://meqsilhouette.readthedocs.io/en/focalpy38/>`_. For the older Python 2 version, `click here <https://meqsilhouette.readthedocs.io/en/v2.7.1/>`_.
+The easiest way to run MeqSilhouette (focalpy38 -- Ubuntu 20.04 + Python 3.8) is to pull the docker image from Docker Hub.
 
-Ubuntu 18.04 + Python 2.7
--------------------------
+Once you have `Docker <https://www.docker.com/>`_ installed on your Ubuntu system, run::
+
+   $ docker pull iniyannatarajan/meqsilhouette:focalpy38
 
-It is recommended to install the dependencies via the `KERN-6 <https://kernsuite.info>`_ software suite::
+Instructions for building a new Docker image from the Dockerfile provided with the source can be found in the relevant section below.
+If you do not have/want Docker, try one of the other options below.
+
+Ubuntu 20.04 + Python 3.8
+-------------------------
+
+It is recommended to install the dependencies via the `KERN-7 <https://kernsuite.info>`_ software suite::
 
    $ sudo apt-get install software-properties-common
-   $ sudo add-apt-repository -s ppa:kernsuite/kern-6
+   $ sudo add-apt-repository -s ppa:kernsuite/kern-7
    $ sudo apt-add-repository multiverse
    $ sudo apt-add-repository restricted
    $ sudo apt-get update
 
 Install the following dependencies via *apt-get*::
 
-   $ sudo apt-get install meqtrees meqtrees-timba tigger tigger-lsm python-astro-tigger \
-   python-astro-tigger-lsm casalite wsclean pyxis python-casacore
+   $ sudo apt-get install meqtrees meqtrees-timba tigger-lsm python3-astro-tigger \
+   python3-astro-tigger-lsm casalite wsclean pyxis python3-casacore
 
 .. note:: The casacore data must be kept up-to-date. This can be done by following the instructions on the `CASA website <https://casaguides.nrao.edu/index.php/Fixing_out_of_date_TAI_UTC_tables_(missing_information_on_leap_seconds)>`_.
 
 Optionally, install Latex (for creating paper-quality plots)::
 
-  $ sudo apt-get install texlive-latex-extra texlive-fonts-recommended dvipng
+  $ sudo apt-get install texlive-latex-extra texlive-fonts-recommended dvipng cm-super
 
-*AATM v0.5* can be obtained from `here <http://www.mrao.cam.ac.uk/~bn204/soft/aatm-0.5.tar.gz>`_. AATM cannot create the executables necessary for running MeqSilhouette without the *boost* libraries. In Ubuntu 18.04, ensure that the packages *libboost-program-options-dev*, *libboost-program-options1.65-dev*, and *libboost-program-options1.65.1* are installed. Once these are installed, proceed as follows::
+*AATM v0.5* can be obtained from `here <http://www.mrao.cam.ac.uk/~bn204/soft/aatm-0.5.tar.gz>`_. AATM cannot create the executables necessary for running MeqSilhouette without the *boost* libraries. In Ubuntu 20.04 install *libboost-program-options-dev* using apt-get. Once this is installed, proceed as follows::
 
    $ cd /path/to/aatm-source-code
    $ ./configure --prefix=/path/to/aatm-installation
@@ -70,15 +77,7 @@ information, simply rebuild the image to eliminate this warning thrown by *CASA*
 Building Docker image
 ---------------------
 
-*Docker* is also supported. Docker can be installed on your system via *apt-get*. 
-
-Docker images for MeqSilhouette v2 are available on Docker Hub::
-
-    $ docker pull iniyannatarajan/meqsilhouette:<tagname>
-
-where *tagname* is the required tag/version name from Docker Hub.
-
-Alternately, build the docker image as follows::
+*Docker* is also supported. Docker can be installed on your system via *apt-get* and the docker image can be built as follows::
 
    $ cd /path/to/Dockerfile
    $ docker build -t meqsilhouette .
@@ -88,6 +87,8 @@ As before, the build process ensures that *casacore* data are up-to-date.
 Known installation issues
 -------------------------
 
+.. note:: One or more of the following issues were originally encountered with Ubuntu 18.04 + Python 2.7 and may possibly be never encountered in the future.
+
 1. If MeqTrees cannot see the *TiggerSkyModel* module that ought to load when *turbo-sim.py* is run (i.e. when an ASCII sky model is used), the parent directory of *Tigger* must be added to PYTHONPATH. Bear in mind that this may cause python version conflicts with other packages. In that case, it is recommended to have Tigger installed in a separate directory such as /opt/Tigger. For manual installation of `Tigger <https://github.com/ska-sa/tigger>`_ and `tigger-lsm <https://github.com/ska-sa/tigger-lsm>`_, refer to their respective repositories. Without this, MeqSilhouette will still work with FITS images as input sky models.
 
 2. If MeqSilhouette cannot find aatm, modify LD_LIBRARY_PATH as follows::

meqsilhouette/__init__.py

@@ -1 +1 @@
-__version__ = "2.7"
+__version__ = "3.0"

meqsilhouette/data/eht230.json

@@ -1,7 +1,7 @@
 {
-"outdirname":"/meqsdata/output/ehttest",
-"input_fitsimage":"/meqsdata/sky_models/singlept",
-"input_fitspol":1,
+"outdirname":"/meqsdata/output/wsclean_timevar_point_allcorrupt",
+"input_fitsimage":"/meqsdata/sky_models/timevar_point",
+"input_fitspol":0,
 "input_changroups":1,
 "output_to_logfile":0,
 "add_thermal_noise":1,

meqsilhouette/driver/run_meqsilhouette.py

@@ -12,7 +12,7 @@
 import pyrap.tables as pt
 import meqsilhouette.framework
 from meqsilhouette.framework.process_input_config import setup_keyword_dictionary, load_json_parameters_into_dictionary
-from meqsilhouette.framework.create_ms import create_ms
+from meqsilhouette.framework.create_ms import create_ms, create_msv2
 from meqsilhouette.framework.SimCoordinator import SimCoordinator
 from meqsilhouette.framework.meqtrees_funcs import make_dirty_image_lwimager
 from meqsilhouette.framework.comm_functions import *
@@ -162,11 +162,12 @@ def run_meqsilhouette(config=None):
     elif parameters['corr_quantbits'] == 2: corr_eff = 0.88
     else: abort('Invalid number of bits used for quantization. Value of "corr_quantbits" in input json file must be 1 or 2')
 
-    info('Creating empty MS with simms')
-    create_ms(MS, input_fitsimage, ms_dict)
+    info('Creating empty MS')
+    #create_ms(MS, input_fitsimage, ms_dict)
+    create_msv2(MS, input_fitsimage, ms_dict)
 
     # Move simms log into the output directory
-    os.system('mv %s %s'%('log-simms.txt', input_copy_path.rsplit('/',1)[0]))
+    #os.system('mv %s %s'%('log-simms.txt', input_copy_path.rsplit('/',1)[0]))
 
     # INI: Write mount types into the MOUNT column in the empty MS prior to generating synthetic data.
     station_mount_types = np.loadtxt(parameters['station_info'], usecols=[19], dtype=str, skiprows=1)
@@ -234,8 +235,9 @@ def run_meqsilhouette(config=None):
             sim_coord.trop_calc_fixdelay_phase_offsets()
             combined_phase_errors += sim_coord.fixdelay_phase_errors
 
-        info('TROPOSPHERE: applying desired combination of phase errors...')
-        sim_coord.apply_phase_errors(combined_phase_errors) 
+        if not isinstance(combined_phase_errors, int):
+            info('TROPOSPHERE: applying desired combination of phase errors...')
+            sim_coord.apply_phase_errors(combined_phase_errors) 
 
         info('All selected tropospheric corruptions applied.')
 

meqsilhouette/framework/SimCoordinator.py

@@ -102,14 +102,13 @@ def __init__(self, msname, output_column, input_fitsimage, input_fitspol, input_
         ### INI: Oversampling factor to use for visibility prediction
         self.oversampling = predict_oversampling
         if predict_seed != -1:
-            self.rng_predict = np.random.RandomState(predict_seed)
+            self.rng_predict = np.random.default_rng(predict_seed)
         else:
-            self.rng_predict = np.random.RandomState()
-
+            self.rng_predict = np.random.default_rng()
         if atm_seed != -1:
-            self.rng_atm = np.random.RandomState(atm_seed)
+            self.rng_atm = np.random.default_rng(atm_seed)
         else:
-            self.rng_atm = np.random.RandomState()
+            self.rng_atm = np.random.default_rng()
 
         ### INI: populate WEIGHT and SIGMA columns
         self.thermal_noise_enabled = thermal_noise_enabled
@@ -197,7 +196,7 @@ def interferometric_sim(self):
             startvis = 0
             endvis = self.vis_per_image
             # INI: cycle through images and simulate including polarisation info, if present.
-            for img_ind in range(self.num_images):
+            for img_ind in range(int(self.num_images)):
                 temp_input_fits = '%s/t%04d'%(self.input_fitsimage,img_ind)
                 info('Simulating visibilities (corr dumps) from %d to %d using input sky model %s'%(startvis,endvis,temp_input_fits))
                 run_wsclean(temp_input_fits, self.input_fitspol, self.input_changroups, startvis, endvis, self.oversampling)
@@ -431,9 +430,12 @@ def trop_return_opacity_sky_temp(self):
             #output = subprocess.check_output(self.ATM_absorption_string(self.chan_freq[0]-self.chan_width, self.chan_freq[-1], self.chan_width,self.average_pwv[ant], self.average_gpress[ant],self.average_gtemp[ant]),shell=True)
             output = subprocess.check_output(ATM_abs_string,shell=True)
             atmfile = open(II('$OUTDIR')+'/atm_output/ATMstring_ant%i.txt'%ant,'w')
-            print>>atmfile,ATM_abs_string
+            print(ATM_abs_string, file=atmfile)
             atmfile.close()
-            atm_abs = file(II('$OUTDIR')+'/atm_output/%satm_abs.txt' % ant, 'w'); atm_abs.write(output); atm_abs.close()
+
+            with open(II('$OUTDIR')+'/atm_output/%satm_abs.txt'%ant, 'wb') as atm_abs:
+                atm_abs.write(output)
+
             if self.num_chan == 1:
                 freq_atm, dry, wet, temp_atm = np.swapaxes(np.expand_dims(np.loadtxt(II('$OUTDIR')+'/atm_output/%satm_abs.txt'%ant, skiprows=1, usecols=[0, 1, 2, 3], delimiter=', \t'), axis=0), 0, 1)
             else:
@@ -509,7 +511,7 @@ def trop_generate_turbulence_phase_errors(self):
             L = np.linalg.cholesky(covmatS) # Cholesky factorise the covariance matrix
 
             # INI: generate random walk error term
-            turb_phase_errors[:, 0, ant] = np.sqrt(1/np.sin(self.elevation_tropshape[:, 0, ant])) * L.dot(self.rng_atm.randn(self.time_unique.shape[0]))
+            turb_phase_errors[:, 0, ant] = np.sqrt(1/np.sin(self.elevation_tropshape[:, 0, ant])) * L.dot(self.rng_atm.standard_normal(self.time_unique.shape[0]))
             turb_phase_errors[:, :, ant] = np.multiply(turb_phase_errors[:, 0, ant].reshape((self.time_unique.shape[0], 1)), (self.chan_freq/self.chan_freq[0]).reshape((1, self.chan_freq.shape[0])))
 
         self.turb_phase_errors = turb_phase_errors
@@ -546,9 +548,12 @@ def trop_ATM_dispersion(self):
 
             output = subprocess.check_output(ATM_disp_string,shell=True)
             atmfile = open(II('$OUTDIR')+'/atm_output/ATMstring_ant%i.txt'%ant,'a')
-            print>>atmfile,ATM_disp_string
+            print(ATM_disp_string, file=atmfile)
             atmfile.close()
-            atm_disp = file(II('$OUTDIR')+'/atm_output/%satm_disp.txt' % ant, 'w'); atm_disp.write(output); atm_disp.close()
+
+            with open(II('$OUTDIR')+'/atm_output/%satm_disp.txt'%ant, 'wb') as atm_disp:
+                atm_disp.write(output)
+
             if self.num_chan == 1:
                 wet_non_disp, wet_disp, dry_non_disp = np.swapaxes(np.expand_dims(np.genfromtxt(II('$OUTDIR')+'/atm_output/%satm_disp.txt'%ant, skip_header=1, usecols=[1, 2, 3], delimiter=',',autostrip=True), axis=0), 0, 1)
             else:
@@ -760,7 +765,7 @@ def pointing_constant_offset(self,pointing_rms, pointing_timescale,PB_FWHM230):
             if self.num_mispoint_epochs != len(self.mjd_ptg_epoch_timecentroid):
                 self.mjd_ptg_epoch_timecentroid = self.mjd_ptg_epoch_timecentroid[:-1]
 
-            self.pointing_offsets = pointing_rms.reshape(self.Nant,1) * self.rng_predict.randn(self.Nant,self.num_mispoint_epochs) # units: arcsec
+            self.pointing_offsets = pointing_rms.reshape(self.Nant,1) * self.rng_predict.standard_normal((self.Nant,self.num_mispoint_epochs)) # units: arcsec
             for ant in range(self.Nant):
                 ind = (self.mjd_ptg_epoch_timecentroid < self.mjd_ant_rise[ant]) \
                     | (self.mjd_ptg_epoch_timecentroid > self.mjd_ant_set[ant])
@@ -834,7 +839,7 @@ def plot_pointing_errors(self):
         #color.cycle_cmap(self.Nant, cmap=cmap) # INI: deprecated
         marker = itertools.cycle(('.', 'o', 'v', '^', 's', '+', '*', 'h', 'D'))
         for i in range(self.Nant):
-            pl.plot(abs(self.pointing_offsets[i,:]),self.pointing_amp_errors[i,:], marker=marker.next(), linestyle='', alpha=1,label=self.station_names[i])
+            pl.plot(abs(self.pointing_offsets[i,:]),self.pointing_amp_errors[i,:], marker=next(marker), linestyle='', alpha=1,label=self.station_names[i])
         pl.xlim(0,np.nanmax(abs(self.pointing_offsets))*1.1)
         pl.ylim(np.nanmin(self.pointing_amp_errors[i,:])*0.8,1.04)
         pl.xlabel(r'Pointing offset, $\rho$ / arcsec', fontsize=FSIZE)
@@ -879,10 +884,10 @@ def add_bjones_manual(self):
             #bjones_interpolated[ant] = spl(self.chan_freq)
             temp_amplitudes_r = spl_r(self.chan_freq)
             temp_amplitudes_l = spl_l(self.chan_freq)
-            temp_phases_r = np.deg2rad(60*self.rng_predict.random_sample(temp_amplitudes_r.shape[0]) - 30) # add random phases between -30 deg to +30 deg
-            temp_phases_l = np.deg2rad(60*self.rng_predict.random_sample(temp_amplitudes_l.shape[0]) - 30) # add random phases between -30 deg to +30 deg
-            self.bjones_interpolated[ant,:,0,0] = np.array(map(cmath.rect, temp_amplitudes_r, temp_phases_r))
-            self.bjones_interpolated[ant,:,1,1] = np.array(map(cmath.rect, temp_amplitudes_l, temp_phases_l))
+            temp_phases_r = np.deg2rad(60*self.rng_predict.random(temp_amplitudes_r.shape[0]) - 30) # add random phases between -30 deg to +30 deg
+            temp_phases_l = np.deg2rad(60*self.rng_predict.random(temp_amplitudes_l.shape[0]) - 30) # add random phases between -30 deg to +30 deg
+            self.bjones_interpolated[ant,:,0,0] = np.array(list(map(cmath.rect, temp_amplitudes_r, temp_phases_r)))
+            self.bjones_interpolated[ant,:,1,1] = np.array(list(map(cmath.rect, temp_amplitudes_l, temp_phases_l)))
 
         # INI: Write the bandpass gains
         np.save(II('$OUTDIR')+'/bterms_timestamp_%d'%(self.timestamp), self.bjones_interpolated)
@@ -1099,7 +1104,7 @@ def make_ms_plots(self):
         #from mpltools import color
         #cmap = pl.cm.Set1
         #color.cycle_cmap(self.Nant, cmap=cmap) # INI: deprecated; use prop_cycle
-        colors = [pl.cm.Set1(i) for i in np.linspace(0, 1, self.nbl)]
+        colors = [pl.cm.Set1(i) for i in np.linspace(0, 1, int(self.nbl))]
         fig, ax = pl.subplots()
         ax.set_prop_cycle(cycler('color', colors))
         for ant0 in range(self.Nant):
@@ -1331,8 +1336,8 @@ def uas2uvdist(ang):
         ax1 = fig.add_subplot(111)
         ax2 = ax1.twiny()
         #x_vlba,y_vlba = np.loadtxt('/home/deane/git-repos/vlbi-sim/output/XMM-LSS/vlba_xmmlss_sigma_vs_uvbin.txt').T #/home/deane/git-repos/vlbi-sim/output/VLBA_COSMOS/vlba_sigma_vs_uvbin.txt',comments='#').T
-        x = np.ravel(zip(uvbins_edges[:-1],uvbins_edges[1:]))
-        y = np.ravel(zip(stdbins,stdbins))
+        x = np.ravel(list(zip(uvbins_edges[:-1],uvbins_edges[1:])))
+        y = np.ravel(list(zip(stdbins,stdbins)))
         #y_minus1ant = np.ravel(zip(stdbins_minus1ant,stdbins_minus1ant))
 
         #ax1.plot(x_vlba,y_vlba*1e6,color='grey',alpha=1,label='VLBA',lw=3)