diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs index e706d6ae8230e08fc4fa8b86e013ee8f63fa4568..beabefc374068ce2d09b42ba4af1c8fd07d49171 100644 --- a/.git-blame-ignore-revs +++ b/.git-blame-ignore-revs @@ -4,4 +4,6 @@ # Migrate code style to Black 70b0c3a8d096f33958f938333a450f186453497f -52418482d3a6906acbf136b9016d157f2d6597ee \ No newline at end of file +52418482d3a6906acbf136b9016d157f2d6597ee +# Migrate using python-lint from standard Makefiles +0ec746a91586bff09d14dc50d082d39d16ae5f0d \ No newline at end of file diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index b5a0d26e81f517608c273cebc60b079ab3f07a4c..bf397c4e69bab38915d73d9f525c070214cdc0a5 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -1,76 +1,37 @@ +image: $SKA_K8S_TOOLS_DOCKER_BUILDER_IMAGE + +variables: + GIT_SUBMODULE_STRATEGY: recursive + stages: + - lint + - build - test + - pages -linting: - stage: test - image: artefact.skao.int/ska-tango-images-pytango-builder:9.3.11 - when: always - before_script: - - pip3 install --upgrade pip - - pip install -r requirements-test.txt - - pip install black - script: - - make lint - after_script: - - mkdir -p ./build/reports - - mv linting.xml ./build/reports - artifacts: - paths: - - build - expire_in: 1 week - when: always - -test: - stage: test - image: python:3.9 - when: always - before_script: - - apt-get update - - apt-get -y install gfortran libboost-dev libboost-python-dev libboost-numpy-dev - - pip install --upgrade pip - - pip install -r requirements.txt -r requirements-test.txt - - export PATH=/home/rascil/.local/bin:$PATH - - make get_rascil_lfs - - cp -r rascil_data/data /usr/local/lib/python3.9/site-packages/ - script: - - make test_unittest - after_script: - - mkdir -p ./build/reports - - mv /builds/ska-telescope/sim/ska-sim-low/*.xml build/reports - coverage: '/^TOTAL.+?(\d+\%)$/' - artifacts: - paths: - - build - expire_in: 1 week - when: always +include: + - project: 'ska-telescope/templates-repository' + file: 'gitlab-ci/includes/python-lint.gitlab-ci.yml' + - project: 'ska-telescope/templates-repository' + file: 'gitlab-ci/includes/python-test.gitlab-ci.yml' + - project: 'ska-telescope/templates-repository' + file: 'gitlab-ci/includes/docs.gitlab-ci.yml' -# "pages" is a special job that is used to upload static content to GitLab -# that can be used to serve a website (i.e. docs for ska-sim-low) (taken from -# https://gitlab.com/ska-telescope/developer.skatelescope.org/-/blob/master/.gitlab-ci.yml) -pages: - stage: test - image: python:3.9 - script: - - apt-get update - - apt-get -y install gfortran libboost-dev libboost-python-dev libboost-numpy-dev - - pip3 install --upgrade pip - - pip install -r requirements.txt -r requirements-test.txt -r docs/requirements-docs.txt - - make get_rascil_lfs - - cp -r rascil_data/data /usr/local/lib/python3.9/site-packages/ - - make -C docs/src html - # for some reason, after the above command, make leaves the working directory - # so in the next cp command we need to specify the full path - - mkdir -p public - - cp -r /builds/ska-telescope/sim/ska-sim-low/docs/build/html/* public - artifacts: - paths: # for this stage to work with "pages", data have to be in "public" - - public - expire_in: 1 week + - project: 'ska-telescope/templates-repository' + file: 'gitlab-ci/includes/finaliser.gitlab-ci.yml' +python-lint: + before_script: + - apt-get update && apt-get -y install cmake + - pip3 install --upgrade pip && pip install -r requirements-test.txt + rules: + - when: always -# Create Gitlab CI badges from CI metrics -# https://developer.skatelescope.org/en/latest/tools/continuousintegration.html#automated-collection-of-ci-health-metrics-as-part-of-the-ci-pipeline -include: - - project: 'ska-telescope/templates-repository' - file: 'gitlab-ci/includes/post_step.yml' +python-test: + image: registry.gitlab.com/ska-telescope/external/rascil/rascil-full:latest + before_script: + - apt-get update && apt-get -y install cmake + - pip3 install --upgrade pip && pip install -r requirements.txt -r requirements-test.txt + rules: + - when: always \ No newline at end of file diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000000000000000000000000000000000000..87158c67c02ef06c88e780a8098e263096fb57eb --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule ".make"] + path = .make + url = https://gitlab.com/ska-telescope/sdi/ska-cicd-makefile.git diff --git a/.make b/.make new file mode 160000 index 0000000000000000000000000000000000000000..a4674571b9b2a0658194c83154bc6b1731b0f144 --- /dev/null +++ b/.make @@ -0,0 +1 @@ +Subproject commit a4674571b9b2a0658194c83154bc6b1731b0f144 diff --git a/.readthedocs.yml b/.readthedocs.yml index dae4495bf54ef39cce593f41796de24f946800ef..f5c2785e904c8692067e06d14cfccf3e4cae32fd 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -24,5 +24,5 @@ build: python: version: 3.7 install: - - requirements: docs/requirements-docs.txt + - requirements: docs/requirements.txt - requirements: requirements.txt \ No newline at end of file diff --git a/Makefile b/Makefile index dc76b62a449a9c77f5f492bf0e9f720f47e06127..6d23892a895f8b05cf71d68d09184e1a12f5b256 100644 --- a/Makefile +++ b/Makefile @@ -1,28 +1,19 @@ -#WORKDIR=./test_results +include .make/base.mk +include .make/python.mk +include .make/docs.mk -#IMAGES=nexus.engageska-portugal.pt/rascil-docker +PYTHON_LINT_TARGET = rfi/ tests/ +PYTHON_SWITCHES_FOR_BLACK = --line-length 79 scripts/ +# we are ignoring all of the pylint errors +PYTHON_SWITCHES_FOR_PYLINT = --exit-zero -lint: -# outputs linting.xml - pylint --exit-zero --output-format=pylint2junit.JunitReporter rfi > linting.xml - pylint --exit-zero --output-format=parseable rfi - black --check . +# E203: whitespace before operations character +# E402: module level import not at top of file (because of matplotlib.use('Agg')) +# W503: line break before binary operator +PYTHON_SWITCHES_FOR_FLAKE8 = --ignore=E203,E402,W503 -test_unittest: - py.test tests \ - --verbose \ - --junitxml=unit-tests.xml \ - --cov=rfi \ - --cov-report=term \ - --cov-report=html:coverage \ - --cov-report=xml:code-coverage.xml \ - --cov-append \ - --cov-branch - -#everything: rm_all build_all test_all - -#.phony: build_all rm_all test_all test_all test_base test_full test_notebook test_all_singularity \ -# test_all_singularity test_base_singularity test_full_singularity test_notebook_singularity everything +docs-pre-build: + pip install -r requirements.txt docs: ## build docs # Outputs docs/build/html diff --git a/README.md b/README.md index 7581e71e439cc651a0407e66c9b391c3c831ef91..cce40ad9eaff002e47d3e5fc296ab79c60d59b82 100644 --- a/README.md +++ b/README.md @@ -7,8 +7,20 @@ Collection of scripts used for various SKA-Low simulations. Please refer to the ## Contribute to this repository -We use [Black](https://github.com/psf/black) to keep the python code style in good shape. -Please make sure you black-formatted your code before merging to master. +[Black](https://github.com/psf/black), [isort](https://pycqa.github.io/isort/), +and various linting tools are used to keep the Python code in good shape. +Please check that your code follows the formatting rules before committing it +to the repository. You can apply Black and isort to the code with: -The linting step in the CI pipeline checks that the code complies with black formatting style, -and fails if that is not the case. \ No newline at end of file +```bash +make python-format +``` + +and you can run the linting checks locally using: + +```bash +make python-lint +``` + +The linting job in the CI pipeline does the same checks, and it will fail if +the code does not pass all of them. Note: we currently ignore pylint failures. diff --git a/docs/requirements-docs.txt b/docs/requirements.txt similarity index 100% rename from docs/requirements-docs.txt rename to docs/requirements.txt diff --git a/docs/src/conf.py b/docs/src/conf.py index b3e934264de775631ddb99110e7f42d0556a2c1b..ca753647e5fb466da954b319ae4170532cbec8f1 100644 --- a/docs/src/conf.py +++ b/docs/src/conf.py @@ -19,7 +19,6 @@ # import os import sys -from unittest import mock sys.path.insert(0, os.path.abspath("../..")) diff --git a/docs/src/lowlevel-rfi/cli.rst b/docs/src/lowlevel-rfi/cli.rst index 81d85651c7bef4eee73522a0999b5d1c583808ed..f5507999c7e7e02a49742d8fd32b3d876b6bce95 100644 --- a/docs/src/lowlevel-rfi/cli.rst +++ b/docs/src/lowlevel-rfi/cli.rst @@ -12,7 +12,7 @@ Pycraf Python script -------------------- .. argparse:: - :filename: ../../rfi/pycraf_scripts/SKA_low_RFI_propagation.py + :filename: rfi/pycraf_scripts/SKA_low_RFI_propagation.py :func: cli_parser :prog: SKA_low_RFI_propagation.py :nodefault: @@ -23,7 +23,7 @@ OSKAR Python script ------------------- .. argparse:: - :filename: ../../rfi/oskar_sim_beam_gain_sing.py + :filename: rfi/oskar_sim_beam_gain_sing.py :func: cli_parser :prog: oskar_sim_beam_gain_sing.py :nodefault: @@ -34,7 +34,7 @@ RASCIL Python script -------------------- .. argparse:: - :filename: ../../rfi/rascil_scripts/simulate_low_rfi_visibility_propagation.py + :filename: rfi/rascil_scripts/simulate_low_rfi_visibility_propagation.py :func: cli_parser :prog: simulate_low_rfi_visibility_propagation.py :nodefault: @@ -45,7 +45,7 @@ Power Spectrum Python script ---------------------------- .. argparse:: - :filename: ../../rfi/rascil_scripts/power_spectrum.py + :filename: rfi/rascil_scripts/power_spectrum.py :func: cli_parser :prog: power_spectrum.py :nodefault: diff --git a/requirements-test.txt b/requirements-test.txt index 488b1b80fa512b13054bf8254310b328ca6ac9ab..d03a9b498d2406d8cffbe7a44c48f3bc935ce34d 100644 --- a/requirements-test.txt +++ b/requirements-test.txt @@ -1,4 +1,6 @@ pylint -pylint2junit +pylint-junit pytest -pytest-cov \ No newline at end of file +pytest-cov +black +isort \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index dc941d455e71c63e00472f8d04fcb7ba926d0269..c710d2cd1594da4009ec77be27f7bb23313daf78 100644 --- a/requirements.txt +++ b/requirements.txt @@ -19,6 +19,3 @@ docopt h5py git+https://gitlab.com/ska-telescope/external/rascil.git@master#egg=rascil - -# pycraf needs pytest to be installed to work... -pytest \ No newline at end of file diff --git a/rfi/common_tools.py b/rfi/common_tools.py index 95d163f02a66bcd0b2ea8f0de27549602591a25f..3bac8c7c3f3a23936a5e4509af59d6e5beb86901 100644 --- a/rfi/common_tools.py +++ b/rfi/common_tools.py @@ -1,9 +1,8 @@ import logging import os - -import numpy as np from typing import Tuple +import numpy as np from astropy import units as u from astropy.coordinates import EarthLocation from rascil.processing_components import ( @@ -14,13 +13,16 @@ from rascil.processing_components import ( LOG = logging.getLogger("rfi-logger") -def default_or_antenna_value(property_name, set_value, default_value, antenna_data): +def default_or_antenna_value( + property_name, set_value, default_value, antenna_data +): """ Return either the user defined value or antenna-based value of a property :param property_name: property name to get value for. Currently either "bandwidth" or "central frequency" - :param set_value: if True, use default value, else use one obtained from antenna_data + :param set_value: if True, use default value, + else use one obtained from antenna_data :param default_value: value from default arguments :param antenna_data: DTV antenna information (pandas.Series) @@ -36,7 +38,8 @@ def default_or_antenna_value(property_name, set_value, default_value, antenna_da except KeyError: raise KeyError( f"The property you provided ({property_name}) is not valid. " - f"Currently only the following can be used: {column_property_dict.keys()}" + f"Currently only the following can be used: " + f"{column_property_dict.keys()}" ) if set_value == "True": @@ -46,7 +49,9 @@ def default_or_antenna_value(property_name, set_value, default_value, antenna_da else: try: value = antenna_data[column_property_dict[property_name]] - LOG.info("%s from csv is %.2f MHz", property_name.capitalize(), value) + LOG.info( + "%s from csv is %.2f MHz", property_name.capitalize(), value + ) except KeyError: value = default_value @@ -65,11 +70,12 @@ def calculate_bandwidth_channel_range(bandwidth, frequency, tx_freq): based on the given bandwidth :param bandwidth: bandwidth to perform calculations for in MHz - :param frequency: numpy array of frequencies between original start and end value + :param frequency: numpy array of frequencies between + original start and end value :param tx_freq: central frequency of the transmitter - :return: the new frequency range values: number of channels, start and end frequencies, - and the array of new frequencies + :return: the new frequency range values: number of channels, + start and end frequencies, and the array of new frequencies """ freq_start = tx_freq - (bandwidth / 2.0) freq_end = tx_freq + (bandwidth / 2.0) @@ -80,8 +86,9 @@ def calculate_bandwidth_channel_range(bandwidth, frequency, tx_freq): freq_start = new_freq[0] except IndexError: LOG.warning( - "Start and end frequencies calculated from bandwidth and central frequency of transmitter" - "are out of range of full frequency range we are simulating. " + "Start and end frequencies calculated from bandwidth " + "and central frequency of transmitter are out of range " + "of full frequency range we are simulating. " "Cannot obtain bandwidth-truncated frequency range." ) return None, None, None, None @@ -92,11 +99,14 @@ def calculate_bandwidth_channel_range(bandwidth, frequency, tx_freq): return n_freqs, freq_start, freq_end, new_freq -def generate_ska_antenna_configuration(use_antfile, antenna_file, rmax, log=LOG): +def generate_ska_antenna_configuration( + use_antfile, antenna_file, rmax, log=LOG +): """ Generate antenna configuration. - :param args: user-provided arguments, must contain: use_antfile, antenna_file + :param args: user-provided arguments, + must contain: use_antfile, antenna_file :param log: logging object also used in Dask :param rmax: maximum distance of station from SKA antenna centre """ diff --git a/rfi/old_run_checker.py b/rfi/old_run_checker.py index cc3671a7ccaddc3ae76fb736eab6b15c8be201eb..30211abd1deb818b17b35c768635cece3ab159eb 100644 --- a/rfi/old_run_checker.py +++ b/rfi/old_run_checker.py @@ -4,10 +4,10 @@ args to file and to check new cmd args against old saved args. So that if a new run has the same settings as the old run, it will not run again unnecessarily. """ +import csv import logging import os import sys -import csv LOG = logging.getLogger("rfi-logger") @@ -64,7 +64,8 @@ def check_for_old_run(args, settings_file): if old_settings == new_settings: LOG.info("New settings match old settings: Exiting.") LOG.info( - "If this is incorrect behaviour, please delete '%s' " "and run again.", + "If this is incorrect behaviour, please delete '%s' " + "and run again.", settings_file, ) sys.exit(0) diff --git a/rfi/oskar_sim_beam_gain_sing.py b/rfi/oskar_sim_beam_gain_sing.py index fb08afbdb20a8b2d819a61a512e4773e3a8bf6c5..cc00448bb6d7f1e6378244a71134657a3bc37e0f 100644 --- a/rfi/oskar_sim_beam_gain_sing.py +++ b/rfi/oskar_sim_beam_gain_sing.py @@ -1,15 +1,19 @@ -import os -import subprocess import argparse import logging -import sys +import os import pprint +import subprocess +import sys + import h5py import numpy as np import pandas as pd +from rfi.common_tools import ( + calculate_bandwidth_channel_range, + default_or_antenna_value, +) from rfi.old_run_checker import check_for_old_run, write_settings -from rfi.common_tools import default_or_antenna_value, calculate_bandwidth_channel_range from rfi.rfi_interface.rfi_data_cube import BeamGainDataCube LOG = logging.getLogger("rfi-logger") @@ -21,11 +25,17 @@ def cli_parser(): """ Defines the command line arguments needed to run the script. """ - parser = argparse.ArgumentParser(description="Calculate beam gain for SKA1-Low") + parser = argparse.ArgumentParser( + description="Calculate beam gain for SKA1-Low" + ) # General inputs - parser.add_argument("--ra", type=float, default=0.0, help="Right Ascension (deg)") - parser.add_argument("--declination", type=float, default=-18.0, help="Declination") + parser.add_argument( + "--ra", type=float, default=0.0, help="Right Ascension (deg)" + ) + parser.add_argument( + "--declination", type=float, default=-18.0, help="Declination" + ) parser.add_argument( "--indir", type=str, @@ -44,7 +54,8 @@ def cli_parser(): parser.add_argument( "--telescope_path", type=str, - default="./data/telescope_files/SKA1-LOW_SKO-0000422_Rev3_38m_SKALA4_spot_frequencies.tm", + default="./data/telescope_files/" + "SKA1-LOW_SKO-0000422_Rev3_38m_SKALA4_spot_frequencies.tm", help="Path to telescope model directory", ) @@ -53,22 +64,27 @@ def cli_parser(): "--input_hdf_file", type=str, default="tv_transmitter_attenuation_cube.hdf5", - help="HDF5 file located in --indir, which contains necessary coordinate information for each RFI source." - "If not specified, use individual az/el and transmitter files located in --indir.", + help="HDF5 file located in --indir, which contains " + "necessary coordinate information for each RFI source." + "If not specified, use individual az/el and " + "transmitter files located in --indir.", ) - # Transmitter file and Az/El file inputs + manually set frequency information + # Transmitter file and Az/El file inputs + # and manually set frequency information parser.add_argument( "--transmitters", type=str, default="data/transmitters/Filtered_DTV_list_1_el.csv", - help="CSV file containing transmitter properties; not used with HFD data", + help="CSV file containing transmitter properties; " + "not used with HFD data", ) parser.add_argument( "--set_freq", type=str, default="False", - help="Choose the central frequency with --freq, otherwise read it from the CSV file;" + help="Choose the central frequency with --freq, " + "otherwise read it from the CSV file;" "not used when input is an HDF5 file.", ) parser.add_argument( @@ -81,8 +97,8 @@ def cli_parser(): "--set_bandwidth", type=str, default="False", - help="Choose the bandwidth with --bandwidth, otherwise read it from the CSV file; " - "not used when input is an HDF5 file.", + help="Choose the bandwidth with --bandwidth, otherwise " + "read it from the CSV file; not used when input is an HDF5 file.", ) parser.add_argument( "--bandwidth", @@ -94,7 +110,8 @@ def cli_parser(): "--N_channels", type=int, default=3, - help="Number of frequency channels (must match nchannels_per_chunk for RFI simulation run); " + help="Number of frequency channels (must match " + "nchannels_per_chunk for RFI simulation run); " "not used when input is an HDF5 file.", ) parser.add_argument( @@ -108,28 +125,38 @@ def cli_parser(): "--choose_range", type=str, default="False", - help="use channels over full frequency range given. If False, default to only over specified " - "bandwidth. If full frequency range larger than bandwidth number of output channels will " - "be those within the bandwidth only. Not used when input is an HDF5 file.", + help="use channels over full frequency range given. " + "If False, default to only over specified bandwidth. " + "If full frequency range larger than bandwidth number " + "of output channels will be those within the bandwidth only. " + "Not used when input is an HDF5 file.", ) parser.add_argument( "--beam_gain_out", type=str, default="transmitter_name", - help="Starting name of output beam gain file for each transmitter. Directory and " - "file-type not required. Not used when input is an HDF5 file.", + help="Starting name of output beam gain file " + "for each transmitter. Directory and file-type " + "not required. Not used when input is an HDF5 file.", ) return parser def _generate_default_oskar_settings( - args, dtv_antenna, in_dir, out_dir, telescope_file_path, telescope_file_name + args, + dtv_antenna, + in_dir, + out_dir, + telescope_file_path, + telescope_file_name, ): """ - Generate settings for OSKAR default run. OSKAR run for only the SKA array centre. + Generate settings for OSKAR default run. + OSKAR run for only the SKA array centre. :param args: Result of argparse.ArgumentParser.parse_args() (namespace) - :param dtv_antenna: pandas.Series with antenna information from args.transmitters file + :param dtv_antenna: pandas.Series with antenna information + from args.transmitters file :param in_dir: directory where transmitter Az_El files are stored :param out_dir: directory to store results :param telescope_file_path: path to CSV file of telescope model @@ -151,9 +178,12 @@ def _generate_default_oskar_settings( ) if args.choose_range == "False" and freq_end - freq_start > bandwidth: - n_freqs, freq_start, freq_end, new_freq = calculate_bandwidth_channel_range( - bandwidth, frequency, tx_freq - ) + ( + n_freqs, + freq_start, + freq_end, + new_freq, + ) = calculate_bandwidth_channel_range(bandwidth, frequency, tx_freq) try: freq_inc = new_freq[1] - new_freq[0] except IndexError: @@ -188,17 +218,22 @@ def _generate_default_oskar_settings( "phase_centre_ra_deg": args.ra, "phase_centre_dec_deg": args.declination, "num_time_steps": 1, - "start_time_utc": "01-01-2000 09:32:00.000", # source at ra-dec has to be visible from station at this time + # source at ra-dec has to be visible from station at this time + "start_time_utc": "01-01-2000 09:32:00.000", "length": 1.0, }, "telescope": { - "input_directory": os.path.join(telescope_file_path, telescope_file_name) + "input_directory": os.path.join( + telescope_file_path, telescope_file_name + ) }, "beam_pattern": { "all_stations": True, "coordinate_type": "Sky model", "coordinate_frame": "Horizon", - "sky_model/file": os.path.join(in_dir, dtv_antenna["Name"] + "_AzEl.txt"), + "sky_model/file": os.path.join( + in_dir, dtv_antenna["Name"] + "_AzEl.txt" + ), "root_path": os.path.join(out_dir, out_file), "telescope_outputs/text_file/cross_power_amp": True, }, @@ -223,11 +258,18 @@ def run_default(app, args, in_dir, out_dir, telescope_model, telescope_parent): :param telescope_parent: path to CSV file of telescope model """ # Loop over transmitters. - filtered_transmitter_df = pd.read_csv(args.transmitters, sep=",", low_memory=False) + filtered_transmitter_df = pd.read_csv( + args.transmitters, sep=",", low_memory=False + ) for i, dtv_antenna in filtered_transmitter_df.iterrows(): LOG.info("Transmitter data:\n%s", pprint.pformat(dtv_antenna)) settings = _generate_default_oskar_settings( - args, dtv_antenna, in_dir, out_dir, telescope_parent, telescope_model + args, + dtv_antenna, + in_dir, + out_dir, + telescope_parent, + telescope_model, ) if settings is None: @@ -257,12 +299,14 @@ def _generate_oskar_settings_from_hdf5( telescope_file_name, ): """ - Generate settings for OSKAR from HDF5 file. OSKAR run for each SKA station individually. + Generate settings for OSKAR from HDF5 file. + OSKAR run for each SKA station individually. :param args: Result of argparse.ArgumentParser.parse_args() (namespace) - :param freq_info: tuple of number of channels, start of frequency range, and frequency increment - (n_channels, freq_start, freq_inc) - :param az_el_coords: list of [azimuth, elevation, distance] coordinates for the RFI source + :param freq_info: tuple of number of channels, start of frequency range, + frequency increment (n_channels, freq_start, freq_inc) + :param az_el_coords: list of [azimuth, elevation, distance] + coordinates for the RFI source :param station_index: index of SKA station (int) :param telescope_file_path: path to CSV file of telescope model :param telescope_file_name: CSV file name of telescope model @@ -285,11 +329,14 @@ def _generate_oskar_settings_from_hdf5( "phase_centre_ra_deg": args.ra, "phase_centre_dec_deg": args.declination, "num_time_steps": 1, - "start_time_utc": "01-01-2000 09:32:00.000", # source at ra-dec has to be visible from station at this time + # source at ra-dec has to be visible from station at this time + "start_time_utc": "01-01-2000 09:32:00.000", "length": 1.0, }, "telescope": { - "input_directory": os.path.join(telescope_file_path, telescope_file_name) + "input_directory": os.path.join( + telescope_file_path, telescope_file_name + ) }, "beam_pattern": { "all_stations": False, @@ -297,7 +344,8 @@ def _generate_oskar_settings_from_hdf5( "coordinate_type": "Sky model", "coordinate_frame": "Horizon", "sky_model/file": "tmp_in.txt", - "root_path": f"tmp_out", # output into temp file, which is read back in later during the run + # output into temp file, which is read back in later during the run + "root_path": "tmp_out", "station_outputs/text_file/auto_power": True, }, } @@ -306,10 +354,12 @@ def _generate_oskar_settings_from_hdf5( def _read_oskar_output_data(out_files): """ - OSKAR produces four files per run, each containing beam gain information for every channel. + OSKAR produces four files per run, each containing + beam gain information for every channel. The four files are one each for Stokes I, Q, U, and V - Here, we read this data back from the files and save it in a dict for later use. + Here, we read this data back from the files and + save it in a dict for later use. :param out_files: list of files produced by OSKAR (path+filename) @@ -329,11 +379,14 @@ def _read_oskar_output_data(out_files): return beam_gain_data -def run_with_hdf5(app, args, in_dir, out_dir, telescope_model, telescope_parent): +def run_with_hdf5( + app, args, in_dir, out_dir, telescope_model, telescope_parent +): """ Generate OSKAR setting using the input RFI data HDF5 file - OSKAR is run for each SKA Low station, and for each TV transmitter - - the output (temporary) OSKAR files are read back in and concatenated into a second HDF5 file. + - the output (temporary) OSKAR files are read back in + and concatenated into a second HDF5 file. - temporary files are removed :param app: available OSKAR binary @@ -361,7 +414,8 @@ def run_with_hdf5(app, args, in_dir, out_dir, telescope_model, telescope_parent) freq_info = (n_freqs, freq_start, freq_inc) LOG.info( - "Running OSKAR for %s RFI sources, %s stations, and %s frequency channels", + "Running OSKAR for %s RFI sources, %s stations, " + "and %s frequency channels", len(source_ids), len(station_ids), n_freqs, @@ -369,7 +423,8 @@ def run_with_hdf5(app, args, in_dir, out_dir, telescope_model, telescope_parent) beam_gain_cube = BeamGainDataCube( args.ra, args.declination, - "01-01-2000 09:32:00.000", # time of observation, same as one in OSKAR settings + # time of observation, same as one in OSKAR settings + "01-01-2000 09:32:00.000", freq_channels, source_ids, len(station_ids), @@ -408,7 +463,7 @@ def run_with_hdf5(app, args, in_dir, out_dir, telescope_model, telescope_parent) if error_code != 0: raise RuntimeError("Failed to run OSKAR") - # read output OSKAR files back in and append beam_gain data to array + # read output OSKAR files back in + append beam_gain data to array out_files = [f for f in os.listdir(".") if f.startswith("tmp_out")] data = _read_oskar_output_data(out_files) @@ -460,7 +515,9 @@ def run_oskar_beam(args): if os.path.exists(args.oskar_path): try: # Command will fail with exception if Singularity is not installed. - subprocess.call(["singularity", "--version"], stdout=subprocess.DEVNULL) + subprocess.call( + ["singularity", "--version"], stdout=subprocess.DEVNULL + ) # Construct the Singularity command to use. app = [ @@ -474,7 +531,10 @@ def run_oskar_beam(args): f"{out_dir}:/out_dir", args.oskar_path, ] + app - LOG.info("Running OSKAR using Singularity container: [%s]", " ".join(app)) + LOG.info( + "Running OSKAR using Singularity container: [%s]", + " ".join(app), + ) # Re-assign directories to those bound inside the container. telescope_parent = "/telescope" @@ -487,16 +547,23 @@ def run_oskar_beam(args): else: # If Singularity image is not found, run natively. LOG.info( - "Singularity image '%s' not found: Running OSKAR natively", args.oskar_path + "Singularity image '%s' not found: Running OSKAR natively", + args.oskar_path, ) - # Generate individual bram_gain output files, load individual az-el text files. + # Generate individual bram_gain output files, + # load individual az-el text files. if args.input_hdf_file == "": - run_default(app, args, in_dir, out_dir, telescope_model, telescope_parent) + run_default( + app, args, in_dir, out_dir, telescope_model, telescope_parent + ) - # Generate beam_gain HDF5 file, load rfi HDF5 file with transmitter coordinates + # Generate beam_gain HDF5 file, + # load rfi HDF5 file with transmitter coordinates else: - run_with_hdf5(app, args, in_dir, out_dir, telescope_model, telescope_parent) + run_with_hdf5( + app, args, in_dir, out_dir, telescope_model, telescope_parent + ) # write args of the call to file write_settings(vars(args), settings_file) diff --git a/rfi/pycraf_scripts/SKA_low_RFI_propagation.py b/rfi/pycraf_scripts/SKA_low_RFI_propagation.py index 0faba0ce4ffd66fd34e2160eef51a0b1deb66506..c4da57264294ac20c761de6b32f93a27e907e0b7 100644 --- a/rfi/pycraf_scripts/SKA_low_RFI_propagation.py +++ b/rfi/pycraf_scripts/SKA_low_RFI_propagation.py @@ -1,33 +1,41 @@ """ -Uses pycraf to calculate the attenuation of the input transmitters at each of the SKA antennas for the required -frequencies. This outputs a numpy array (.npy file) for each transmitter in the input csv file. +Uses pycraf to calculate the attenuation of the +input transmitters at each of the SKA antennas for +the required frequencies. This outputs a numpy array +(.npy file) for each transmitter in the input csv file. """ +import argparse import logging import os import pprint import sys -import argparse + import matplotlib matplotlib.use("Agg") +from typing import Optional, Tuple, Union + import matplotlib.pyplot as plt import numpy as np import pandas as pd - +from astropy import ( + units as u, # must be after pycraf to avoid astropy import conflicts +) from numpy import ndarray from pycraf import pathprof -from astropy import units as u # must be after pycraf to avoid astropy import conflicts -from typing import Optional, Union, Tuple -from rfi.pycraf_scripts.SKA_low_pycraf_propagation import calc_prop_atten +from rfi.common_tools import ( + calculate_bandwidth_channel_range, + default_or_antenna_value, +) from rfi.old_run_checker import check_for_old_run, write_settings -from rfi.common_tools import default_or_antenna_value, calculate_bandwidth_channel_range +from rfi.pycraf_scripts.SKA_low_pycraf_propagation import calc_prop_atten from rfi.rfi_interface.rfi_data_cube import ( - RFISource, - RFISignal, - DataCubePerSource, DataCube, + DataCubePerSource, + RFISignal, + RFISource, ) LOG = logging.getLogger("rfi-logger") @@ -47,7 +55,8 @@ def cli_parser(): "--set_freq", type=str, default="False", - help="Choose the central frequency with --freq, otherwise read it from the csv file", + help="Choose the central frequency with --freq, " + "otherwise read it from the csv file", ) parser.add_argument( "--freq", type=float, default=177.5, help="Central frequency (MHz)" @@ -56,14 +65,18 @@ def cli_parser(): "--set_bandwidth", type=str, default="True", - help="Choose the bandwidth with --bandwidth, otherwise read it from the csv file", + help="Choose the bandwidth with --bandwidth, " + "otherwise read it from the csv file", + ) + parser.add_argument( + "--bandwidth", type=float, default=7.0, help="Bandwidth (MHz)" ) - parser.add_argument("--bandwidth", type=float, default=7.0, help="Bandwidth (MHz)") parser.add_argument( "--n_channels", type=int, default=3, - help="Number of frequency channels. (Must match nchannels_per_chunk for RFI simulation run)", + help="Number of frequency channels. " + "(Must match nchannels_per_chunk for RFI simulation run)", ) parser.add_argument( "--frequency_range", @@ -82,26 +95,29 @@ def cli_parser(): "--trans_out", type=str, default="infile", - help="Name of output transmitter list. File-type not required. If 'infile' will supplement the " - "input name of the input transmitter file. If not full path, it will be written to " + help="Name of output transmitter list. File-type not required. " + "If 'infile' will supplement the input name of the input " + "transmitter file. If not full path, it will be written to " "output_dir directory.", ) parser.add_argument( "--non_below_el", type=str, default="True", - help="If transmitter elevation to array centre < 0 deg, remove from list and do not simulate.", + help="If transmitter elevation to array centre < 0 deg, " + "remove from list and do not simulate.", ) parser.add_argument( "--srtm_directory", type=str, default="../data/SRTM_data", - help="Directory for the SRTM files required by pycraf for terrain information.", + help="Directory for the SRTM files required by pycraf " + "for terrain information.", ) parser.add_argument( "--antenna_file", type=str, - default="../data/telescope_files/SKA1-LOW_SKO-0000422_Rev3_38m_SKALA4_spot_frequencies.tm/layout_wgs84.txt", + default="../data/telescope_files/SKA1-LOW_SKO-0000422_Rev3_38m_SKALA4_spot_frequencies.tm/layout_wgs84.txt", # noqa: E501 help="Location of text files with antenna locations", ) parser.add_argument( @@ -111,7 +127,10 @@ def cli_parser(): help="Maximum distance of station from centre (m)", ) parser.add_argument( - "--station_skip", type=int, default=1, help="Decimate stations by this factor" + "--station_skip", + type=int, + default=1, + help="Decimate stations by this factor", ) parser.add_argument( "--output_dir", @@ -123,19 +142,23 @@ def cli_parser(): "--array_centre", type=float, default=["Low_E4", 116.4525771, -26.60055525], - help="List containing name, latitude (degs), longitude (degs) for the SKA Low array centre", + help="List containing name, latitude (degs), " + "longitude (degs) for the SKA Low array centre", ) parser.add_argument( "--plot_attenuation", type=str, default=False, - help="Output plot of attenuation values for each transmitter at the array centre.", + help="Output plot of attenuation values for each " + "transmitter at the array centre.", ) parser.add_argument( "--n_time_chunks", type=int, default=1, - help="Number of time samples to simulate. (Same as the RASCIL-based part's --nintegrations_per_chunk arg.)", + help="Number of time samples to simulate. " + "(Same as the RASCIL-based part's " + "--nintegrations_per_chunk arg.)", ) parser.add_argument( "--frequency_variable", @@ -179,7 +202,8 @@ def cli_parser(): "--height_rg", type=float, default=2.1, - help="Assumed height of receiver above ground (m). See pycraf documentation.", + help="Assumed height of receiver above ground (m). " + "See pycraf documentation.", ) parser.add_argument( "--diam", @@ -191,14 +215,15 @@ def cli_parser(): "--zones", type=str, default=[pathprof.CLUTTER.UNKNOWN, pathprof.CLUTTER.UNKNOWN], - help="List of clutter types for transmitter and receiver, default is unknown. " - "See pycraf documentation.", + help="List of clutter types for transmitter and " + "receiver, default is unknown. See pycraf documentation.", ) parser.add_argument( "--hprof_step", type=float, default=1000, - help="Distance resolution of the calculated solution. See pycraf documentation.", + help="Distance resolution of the calculated solution. " + "See pycraf documentation.", ) return parser @@ -207,7 +232,8 @@ def transmitter_below_horizon(check_horizon, elevation): """ Perform transmitter-horizon checks. - :param check_horizon: if True, check where transmitter is with respect to the horizon + :param check_horizon: if True, check where transmitter + is with respect to the horizon :param elevation: elevation of the transmitter :return bool: True if transmitter is below horizon @@ -231,20 +257,24 @@ def plot_attenuation(atten_ant, freqs, outputdir, tx_name): # TODO: Danielle says that plt didn't work without fig # being defined and fig wouldn't work with plotting. # This needs investigation. - fig = plt.figure(figsize=(10, 10)) + plt.figure(figsize=(10, 10)) plt.plot(freqs, np.transpose(atten_ant)) plt.xlabel("Frequency (MHz)") plt.ylabel("Total attenuation (dB)") plt.title("Station attenuation") plt.grid() plt.savefig( - outputdir + "Station_attenuation_" + tx_name + ".png", bbox_inches="tight" + outputdir + "Station_attenuation_" + tx_name + ".png", + bbox_inches="tight", ) -def construct_csv_file_path(transmitter_file_name, filtered_transmitters, outputdir): +def construct_csv_file_path( + transmitter_file_name, filtered_transmitters, outputdir +): """ - Generate the CSV file name and path where the transmitter information is saved + Generate the CSV file name and path where the + transmitter information is saved :param transmitter_file_name: provided by user, same as args.trans_out :param filtered_transmitters: file containing input transmitter properties @@ -349,19 +379,22 @@ def calculate_apparent_power( n_chans: number of frequency channels :param n_times: number of time samples - :param frequency_range: frequency range within which transmitter emits (float) [MHz] + :param frequency_range: frequency range within which + transmitter emits (float) [MHz] :param transmitter_power: transmitter's emitted power (float) :param attenuation: numpy array of attenuation values [n_ants, n_chan] :param frequency_variable: whether signal should be frequency-variable :param time_variable: whether signal should be time-variable - :return apparent_power: attenuated, apparent power of the transmitter at SKA stations + :return apparent_power: attenuated, apparent power of + the transmitter at SKA stations """ if frequency_range == 0: frequency_range = 1.0e-6 - # this bit is from rascil.processing_components.simulation.rfi.simulate_DTV_prop + # this bit is from + # rascil.processing_components.simulation.rfi.simulate_DTV_prop amp = np.sqrt(transmitter_power / (2.0 * frequency_range * 1e6)) # [W/Hz] n_chans = attenuation.shape[1] @@ -387,10 +420,14 @@ def calculate_apparent_power( else: transmitter_signal += amp - # this is from rascil.processing_components.simulation.rfi.create_propagators_prop - # Danielle thinks that the attenuation information is per m^2, the related calculations done - # within pycraf --> the units of propagation would be 1/m^2 - propagation = np.sqrt(np.power(10, attenuation / -10.0)) # [n_ants, n_chans] + # this is from + # rascil.processing_components.simulation.rfi.create_propagators_prop + # Danielle thinks that the attenuation information is per m^2, + # the related calculations done within pycraf + # --> the units of propagation would be 1/m^2 + propagation = np.sqrt( + np.power(10, attenuation / -10.0) + ) # [n_ants, n_chans] # [n_times, n_ants, n_chans] apparent_power = ( @@ -407,12 +444,15 @@ def rfi_attenuation(args, save_azel_to_file=True, hdf5_only=False): apparent transmitter power at the give SKA stations. :param args: User defined arguments - :param save_azel_to_file: if True, save azimuth and elevation information into a file for OSKAR to user; + :param save_azel_to_file: if True, save azimuth and elevation information + into a file for OSKAR to user; set to False for tests - :param hdf5_only: if True, do not produce any files, but only return data to be saved in HDF5 + :param hdf5_only: if True, do not produce any files, + but only return data to be saved in HDF5 - :return: data cube containing transmitter information, apparent source coordinates, - and apparent transmitter power data + :return: data cube containing transmitter information, + apparent source coordinates, and + apparent transmitter power data """ LOG.info("Starting LOW propagation calculation") @@ -448,7 +488,9 @@ def rfi_attenuation(args, save_azel_to_file=True, hdf5_only=False): rx_name, rx_lon, rx_lat = args.array_centre # load transmitter csv to select licence_no # - filtered_transm_df = pd.read_csv(filtered_transmitters, sep=",", engine="python") + filtered_transm_df = pd.read_csv( + filtered_transmitters, sep=",", engine="python" + ) filtered_transm_df = filtered_transm_df.set_index("Licence Number") transmitter_dict = {} @@ -456,9 +498,10 @@ def rfi_attenuation(args, save_azel_to_file=True, hdf5_only=False): tv_antenna_data = [] for i, dtv_ant in filtered_transm_df.iterrows(): - # need to set them here again, because when a new frequency range is calculated via - # calculate_bandwidth_channel_range and that returns Nones, these values are overwritten - # we want to start with the original values for each transmitter + # need to set them here again, because when a new frequency + # range is calculated via calculate_bandwidth_channel_range + # and that returns Nones, these values are overwritten we want + # to start with the original values for each transmitter n_freqs = args.n_channels freq_start = args.frequency_range[0] freq_end = args.frequency_range[1] @@ -511,10 +554,15 @@ def rfi_attenuation(args, save_azel_to_file=True, hdf5_only=False): continue # Only caclulate attenuation over specified transmitter bandwidth. - # If full frequency range is larger than bandwidth, the number of output channels will - # be those within the bandwidth only. + # If full frequency range is larger than bandwidth, the number + # of output channels will be those within the bandwidth only. if freq_end - freq_start > bandwidth: - n_freqs, freq_start, freq_end, _ = calculate_bandwidth_channel_range( + ( + n_freqs, + freq_start, + freq_end, + _, + ) = calculate_bandwidth_channel_range( bandwidth, frequency, tx_freq ) if n_freqs is None: @@ -534,7 +582,8 @@ def rfi_attenuation(args, save_azel_to_file=True, hdf5_only=False): else: LOG.info("Calculating attenuation for all %d channels", n_freqs) - # Calculate propagation attenuation using pycraf + return [az, el, dist] for each ska station + # Calculate propagation attenuation using pycraf + # + return [az, el, dist] for each ska station atten_ant, freqs, apparent_coords, ska_ids = calc_prop_atten( tx_freq, args.temperature, @@ -622,7 +671,9 @@ def rfi_attenuation(args, save_azel_to_file=True, hdf5_only=False): if not hdf5_only: if save_azel_to_file: fp = open(outputdir + tx_name + "_AzEl.txt", "w") - print(az_el[0], az_el[1], file=fp) # save azel to file for use in OSKAR + print( + az_el[0], az_el[1], file=fp + ) # save azel to file for use in OSKAR fp.close() if not hdf5_only: @@ -631,7 +682,8 @@ def rfi_attenuation(args, save_azel_to_file=True, hdf5_only=False): ) LOG.info( - "%d of %d transmitters below the horizon and removed " "from the output list", + "%d of %d transmitters below the horizon and removed " + "from the output list", len(trans_below), len(filtered_transm_df.index), ) @@ -653,7 +705,11 @@ def main(): times = [None] * args.n_time_chunks cube = DataCube( - times, freqs, station_ids, rmax=args.rmax, station_skip=args.station_skip + times, + freqs, + station_ids, + rmax=args.rmax, + station_skip=args.station_skip, ) for source_data in rfi_results: source_data.time_samples = np.array(times) diff --git a/rfi/pycraf_scripts/SKA_low_pycraf_propagation.py b/rfi/pycraf_scripts/SKA_low_pycraf_propagation.py index 84ff20455ca095564d3387b2b5a1782fe4fa46d9..9c23cbbcf9ee41e7a33ea7225b9acf1a4236987a 100644 --- a/rfi/pycraf_scripts/SKA_low_pycraf_propagation.py +++ b/rfi/pycraf_scripts/SKA_low_pycraf_propagation.py @@ -2,20 +2,25 @@ # coding: utf-8 """ -Core pycraf calculation of the propagation attenuation of a transmitter to each of the SKA low stations defined in the -input workbook. Returns a 2D float array of frequency and attenuation values. +Core pycraf calculation of the propagation attenuation of +a transmitter to each of the SKA low stations defined in +the input workbook. Returns a 2D float array of frequency +and attenuation values. """ # Original by Federico Di Vruno # Alterations by DMF import logging -import numpy as np -from pycraf import pathprof, antenna, geometry -from astropy import units as u # must be after pycraf to avoid astropy import conflicts import time +import numpy as np +from astropy import ( + units as u, # must be after pycraf to avoid astropy import conflicts +) +from pycraf import antenna, geometry, pathprof from rascil.processing_components import select_configuration + from rfi.common_tools import generate_ska_antenna_configuration # allow download of missing SRTM data @@ -43,7 +48,8 @@ def calc_prop_atten( direct, rmax, station_skip=None, - locations_file="./telescope_files/SKA1-LOW_SKO-0000422_Rev3_38m_SKALA4_spot_frequencies.tm/layout_wgs84.txt", + locations_file="./telescope_files/SKA1-LOW_SKO-0000422_Rev3_38m_" + "SKALA4_spot_frequencies.tm/layout_wgs84.txt", polarization=None, srtm_directory="./data/SRTM_data", ): @@ -75,7 +81,7 @@ def calc_prop_atten( - corresponding frequencies; dim = [Nfreq] - array containing azimuth, elevation, distance for each station; dim = [Nant] (i.e. number of SKA stations) and the station ids used for the attenuation calculations - """ + """ # noqa: E501 pathprof.SrtmConf.set(srtm_dir=srtm_directory) freq = freq * u.MHz @@ -91,8 +97,11 @@ def calc_prop_atten( lon_tx, lat_tx = tx_lon * u.deg, tx_lat * u.deg # For LOW, get list of station locations from the txt file - # then sub-select the ones which are at the right place from the array centre - all_long_x, all_lat_x = np.loadtxt(locations_file, delimiter=",", unpack=True) + # then sub-select the ones which are at the right place + # from the array centre + all_long_x, all_lat_x = np.loadtxt( + locations_file, delimiter=",", unpack=True + ) low = generate_ska_antenna_configuration(True, locations_file, rmax) LOG.info("Maximum station distance from array centre: {} m".format(rmax)) @@ -101,10 +110,15 @@ def calc_prop_atten( # apply station skipping sub_low = select_configuration(low, low.names[::station_skip]) - LOG.info("There are {na:d} stations after decimation".format(na=len(sub_low.names))) + LOG.info( + "There are {na:d} stations after decimation".format( + na=len(sub_low.names) + ) + ) - # TODO: differences in amplitude in the final visibility are about an order of magnitude - # need to investigate the reason, and once satisfied, update code to use NEW instead of OLD + # TODO: differences in amplitude in the final visibility + # are about an order of magnitude need to investigate the reason, + # and once satisfied, update code to use NEW instead of OLD # *NEW*: taking the stations that match by name in the configuration # stations_to_use = [int(x.strip("LOW_")) - 1 for x in low.names.values] # long_x = all_long_x[stations_to_use] @@ -115,20 +129,22 @@ def calc_prop_atten( lat_x = all_lat_x[: len(sub_low.names)] n_ant = len(long_x) - # when the coordinates file contains fewer entries than stations in the configuration - # we need to know which of those stations (and how many) could be used + # when the coordinates file contains fewer entries + # than stations in the configuration we need to know + # which of those stations (and how many) could be used # TODO: this will probably be different once we use the proper stations station_ids = sub_low.names.data[:n_ant] - """ - =========================================== - Find the attenuation at each position of SKA antennas using path info for each freq. - =========================================== - """ + + # Find the attenuation at each position of SKA antennas + # using path info for each freq: attenpath_calc_start = time.time() LOG.info( - "Calculating frequency dependent attenuation for " "each station from path..." + "Calculating frequency dependent attenuation for " + "each station from path..." + ) + freqs = ( + np.logspace(np.log10(freq_start), np.log10(freq_end), N_freqs) * u.MHz ) - freqs = np.logspace(np.log10(freq_start), np.log10(freq_end), N_freqs) * u.MHz diameter_fl = diam * u.m @@ -178,7 +194,9 @@ def calc_prop_atten( if direct == "OD": G_eff = G_max_fl else: - G_eff = antenna.fl_pattern(ang_dist, diameter_fl, wavelen_fl, G_max_fl) + G_eff = antenna.fl_pattern( + ang_dist, diameter_fl, wavelen_fl, G_max_fl + ) results = pathprof.atten_path_fast( freqs[i], @@ -194,6 +212,8 @@ def calc_prop_atten( Atten_ant[k, i] = results["L_b"][-1].value - G_eff.value attenpath_calc_end = time.time() - LOG.info("Completed in %f sec." % (attenpath_calc_end - attenpath_calc_start)) + LOG.info( + "Completed in %f sec." % (attenpath_calc_end - attenpath_calc_start) + ) return Atten_ant, freqs, angular_dist_ant, station_ids diff --git a/rfi/pycraf_scripts/supporting_scripts/Test_pycraf_SKA_LOW_antenna.py b/rfi/pycraf_scripts/supporting_scripts/Test_pycraf_SKA_LOW_antenna.py index d49a2c527966bf9650510053f80e7e577cabba5c..3e31bb7753a2cbdd8431dfa7e2ac4f08db2cdb4c 100644 --- a/rfi/pycraf_scripts/supporting_scripts/Test_pycraf_SKA_LOW_antenna.py +++ b/rfi/pycraf_scripts/supporting_scripts/Test_pycraf_SKA_LOW_antenna.py @@ -1,5 +1,6 @@ #!/usr/bin/env python # coding: utf-8 +# flake8: noqa """ Test pycraf using SKA LOW coordinates. @@ -13,13 +14,14 @@ function of frequency. This script is designed to be an test case providing an i Pycraf documentation. """ -from collections import namedtuple, OrderedDict -import numpy as np +import time +from collections import OrderedDict, namedtuple + import matplotlib.pyplot as plt -from pycraf import pathprof -from pycraf import conversions as cnv +import numpy as np from astropy import units as u -import time +from pycraf import conversions as cnv +from pycraf import pathprof # allow download of missing SRTM data: srtm_directory = "../data/SRTM_data/" @@ -100,18 +102,20 @@ Map size inputs ======================================= """ -do_map_solution = True # set to variable if want map solutions to be calculated and map +do_map_solution = ( + True # set to variable if want map solutions to be calculated and map +) # plots produced. Unless set below will provide full resolution plots localised to receivers. doplotAll = True # True for plots and calculated maps to be same area i.e. cover area over transmitter and receiver. -choose_resolution = ( - True # if doplotAll True, by default a full resolution full spatial extent plot -) +choose_resolution = True # if doplotAll True, by default a full resolution full spatial extent plot # (choose_resolution=False) or specify lower resolution options for quicker plotting # (choose_resolution=True), specify options below. -hprof_step = 10000 * u.m # Sets distance resolution for all attenuation calculations. +hprof_step = ( + 10000 * u.m +) # Sets distance resolution for all attenuation calculations. if do_map_solution: # Produce plots locally to receiver at full resolution @@ -162,7 +166,9 @@ if do_map_solution: atten_map_size_lon = map_size_x * u.deg atten_map_size_lat = map_size_y * u.deg atten_map_resolution = 0.1 * u.deg - atten_hprof_step = 15000 * u.m # overrides map_resolution when given + atten_hprof_step = ( + 15000 * u.m + ) # overrides map_resolution when given print("Full distance plot at restricted resolution") else: print("Full distance plot") @@ -181,7 +187,10 @@ omega = 0.0 * u.percent # fraction of path over sea temperature = 290.0 * u.K pressure = 1013.0 * u.hPa timepercent = 0.02 * u.percent # see P.452 for explanation -h_tg, h_rg = 3 * u.m, 2.1 * u.m # height of the receiver and transmitter above gnd +h_tg, h_rg = ( + 3 * u.m, + 2.1 * u.m, +) # height of the receiver and transmitter above gnd G_t, G_r = 0 * cnv.dBi, 0 * cnv.dBi zone_t, zone_r = pathprof.CLUTTER.UNKNOWN, pathprof.CLUTTER.UNKNOWN polarization = 0 # 0 horizontal, 1 vertical @@ -243,7 +252,9 @@ if do_map_solution: vmin, vmax = np.min(_heightmap), np.max(_heightmap) # Produces terrain colormap and norm for plotting - terrain_cmap, terrain_norm = pathprof.terrain_cmap_factory(sealevel=vmin, vmax=vmax) + terrain_cmap, terrain_norm = pathprof.terrain_cmap_factory( + sealevel=vmin, vmax=vmax + ) if doplotAll: @@ -256,11 +267,16 @@ if do_map_solution: else: # Limit to receiver area ax.set_aspect( - abs(map_extent[2] - map_extent[0]) / abs(map_extent[3] - map_extent[1]) + abs(map_extent[2] - map_extent[0]) + / abs(map_extent[3] - map_extent[1]) ) # nextent = (map_extent[0], map_extent[2], map_extent[1], map_extent[3]) - lon_ind = np.where((_lons > map_extent[0]) & (_lons < map_extent[2]))[0] - lat_ind = np.where((_lats > map_extent[1]) & (_lats < map_extent[3]))[0] + lon_ind = np.where((_lons > map_extent[0]) & (_lons < map_extent[2]))[ + 0 + ] + lat_ind = np.where((_lats > map_extent[1]) & (_lats < map_extent[3]))[ + 0 + ] _lons_fin = _lons[lon_ind] _lats_fin = _lats[lat_ind] ax.set_ylim([_lats_fin[0], _lats_fin[-1]]) @@ -394,11 +410,16 @@ if do_map_solution: else: # Limit to receiver area ax.set_aspect( - abs(map_extent[2] - map_extent[0]) / abs(map_extent[3] - map_extent[1]) + abs(map_extent[2] - map_extent[0]) + / abs(map_extent[3] - map_extent[1]) ) # nextent = (map_extent[0], map_extent[2], map_extent[1], map_extent[3]) - lon_ind = np.where((_lons > map_extent[0]) & (_lons < map_extent[2]))[0] - lat_ind = np.where((_lats > map_extent[1]) & (_lats < map_extent[3]))[0] + lon_ind = np.where((_lons > map_extent[0]) & (_lons < map_extent[2]))[ + 0 + ] + lat_ind = np.where((_lats > map_extent[1]) & (_lats < map_extent[3]))[ + 0 + ] _lons_fin = _lons[lon_ind] _lats_fin = _lats[lat_ind] nextent = (_lons_fin[0], _lons_fin[-1], _lats_fin[0], _lats_fin[-1]) @@ -417,7 +438,9 @@ if do_map_solution: # colourbar # cbar = fig.colorbar(cim, cax=cbax, orientation='horizontal', ) - cbar = fig.colorbar(cim, orientation="horizontal", fraction=0.046, pad=0.08) + cbar = fig.colorbar( + cim, orientation="horizontal", fraction=0.046, pad=0.08 + ) cbar.set_label(r"Path propagation loss (dB)", color="black") ctics = np.arange(vmin, vmax, 10) @@ -450,7 +473,9 @@ if do_map_solution: ) ax.scatter(site.coord[0], site.coord[1], marker="o", c="b") - for i in range(len(long_x)): # Puts a mark in every place there is an antenna + for i in range( + len(long_x) + ): # Puts a mark in every place there is an antenna ax.scatter(long_x[i], lat_x[i], marker="o", c="w") plt.title("Attenuation map, Freq = " + str(freq)) @@ -590,7 +615,9 @@ if do_map_solution: """ attenpath_calc_start = time.time() -print("Calculating frequency dependent attenuation for each station from path...") +print( + "Calculating frequency dependent attenuation for each station from path..." +) N_ant = len(long_x) - 1 freqs = np.logspace(np.log10(freq_start), np.log10(freq_end), N_freqs) * u.MHz @@ -695,14 +722,28 @@ print( L_bd ) ) -print("L_bs: {0.value:5.2f} {0.unit} - Tropospheric scatter loss".format(L_bs)) -print("L_ba: {0.value:5.2f} {0.unit} - Ducting/layer reflection loss".format(L_ba)) -print("L_b: {0.value:5.2f} {0.unit} - Complete path propagation loss".format(L_b)) +print( + "L_bs: {0.value:5.2f} {0.unit} - Tropospheric scatter loss".format( + L_bs + ) +) +print( + "L_ba: {0.value:5.2f} {0.unit} - Ducting/layer reflection loss".format( + L_ba + ) +) +print( + "L_b: {0.value:5.2f} {0.unit} - Complete path propagation loss".format( + L_b + ) +) print( "L_b_corr: {0.value:5.2f} {0.unit} - As L_b but with clutter correction".format( L_b_corr ) ) print( - "L: {0.value:5.2f} {0.unit} - As L_b_corr but with gain correction".format(L) + "L: {0.value:5.2f} {0.unit} - As L_b_corr but with gain correction".format( + L + ) ) diff --git a/rfi/rascil_scripts/power_spectrum.py b/rfi/rascil_scripts/power_spectrum.py index f05c7876be7c123af0b42d2248cf9155f09b183a..75fe19f6d49f3fa93eb9c7c01dc205c0afd78c4e 100644 --- a/rfi/rascil_scripts/power_spectrum.py +++ b/rfi/rascil_scripts/power_spectrum.py @@ -1,8 +1,8 @@ # coding: utf-8 """ Power spectrum for an image, adapted from Fred Dulwich code """ -import os import argparse +import os import astropy.constants as consts import matplotlib @@ -10,13 +10,12 @@ import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy - +from rascil.data_models.xarray_coordinate_support import griddata_wcs from rascil.processing_components import ( fft_image_to_griddata, import_image_from_fits, show_image, ) -from rascil.data_models.xarray_coordinate_support import griddata_wcs def radial_profile(image, centre=None): @@ -29,7 +28,9 @@ def radial_profile(image, centre=None): def cli_parser(): - parser = argparse.ArgumentParser(description="Display power spectrum of image") + parser = argparse.ArgumentParser( + description="Display power spectrum of image" + ) parser.add_argument("--image", type=str, default=None, help="Image name") parser.add_argument( "--signal_channel", @@ -38,7 +39,10 @@ def cli_parser(): help="Channel containing both signal and noise", ) parser.add_argument( - "--noise_channel", type=int, default=0, help="Channel containing noise only" + "--noise_channel", + type=int, + default=0, + help="Channel containing noise only", ) parser.add_argument( "--resolution", @@ -82,16 +86,18 @@ def power_spectrum(args): imfft = fft_image_to_griddata(im) print(imfft) - omega = numpy.pi * resolution ** 2 / (4 * numpy.log(2.0)) + omega = numpy.pi * resolution**2 / (4 * numpy.log(2.0)) wavelength = consts.c / numpy.average(im.frequency) - kperjy = 1e-26 * wavelength ** 2 / (2 * consts.k_B * omega) + kperjy = 1e-26 * wavelength**2 / (2 * consts.k_B * omega) im_spectrum = imfft.copy() im_spectrum["pixels"].data = kperjy.value * numpy.abs(imfft["pixels"].data) noisy = numpy.max(im_spectrum["pixels"].data[noise_channel, 0]) > 0.0 profile = radial_profile(im_spectrum["pixels"].data[signal_channel, 0]) - noise_profile = radial_profile(im_spectrum["pixels"].data[noise_channel, 0]) + noise_profile = radial_profile( + im_spectrum["pixels"].data[noise_channel, 0] + ) plt.clf() cellsize_uv = numpy.abs(griddata_wcs(imfft).wcs.cdelt[0]) diff --git a/rfi/rascil_scripts/simulate_low_rfi_visibility_propagation.py b/rfi/rascil_scripts/simulate_low_rfi_visibility_propagation.py index 92ec660d1fb347b6a0a79706ffa2eea20d5939c6..c81a2b271536767c5dc389319a6d9bb4ea65c7b0 100644 --- a/rfi/rascil_scripts/simulate_low_rfi_visibility_propagation.py +++ b/rfi/rascil_scripts/simulate_low_rfi_visibility_propagation.py @@ -33,43 +33,42 @@ The simulation is implemented in some functions in RASCIL, and the script simula is available in the SKA Github repository ska-sim-low. Distributed processing is implemented via Dask. The outputs are fits file and plots of the images: on signal channels and on pure noise channels, and for the source of interest and the Southern Celestial Pole. The un-averaged MeasurementSets are also output, one per time chunk. -""" +""" # noqa: E501 +import argparse import logging import os import pprint import sys import time +import astropy.constants as const import h5py +import matplotlib import numpy -import argparse - from astropy import units as u from astropy.coordinates import SkyCoord -import astropy.constants as const - from rascil.data_models.polarisation import PolarisationFrame +from rascil.processing_components.simulation.configurations import ( + select_configuration, +) from rascil.processing_components.simulation.rfi import ( calculate_averaged_correlation, simulate_rfi_block_prop, ) from rascil.processing_components.util.array_functions import average_chunks -from rascil.processing_components.visibility.base import export_blockvisibility_to_ms -from rascil.workflows.rsexecute.visibility.blockvisibility_rsexecute import ( - concatenate_blockvisibility_time_rsexecute, +from rascil.processing_components.visibility.base import ( + create_blockvisibility, + export_blockvisibility_to_ms, ) from rascil.workflows.rsexecute.execution_support.rsexecute import rsexecute -from rascil.processing_components.simulation.configurations import ( - select_configuration, +from rascil.workflows.rsexecute.visibility.blockvisibility_rsexecute import ( + concatenate_blockvisibility_time_rsexecute, ) -from rascil.processing_components.visibility.base import create_blockvisibility from rfi.common_tools import generate_ska_antenna_configuration from rfi.old_run_checker import check_for_old_run, write_settings -import matplotlib - matplotlib.use("Agg") log = logging.getLogger("rfi-logger") @@ -78,8 +77,12 @@ log.addHandler(logging.StreamHandler(sys.stdout)) def cli_parser(): - parser = argparse.ArgumentParser(description="Simulate RFI data with RASCIL") - parser.add_argument("--seed", type=int, default=18051955, help="Random number seed") + parser = argparse.ArgumentParser( + description="Simulate RFI data with RASCIL" + ) + parser.add_argument( + "--seed", type=int, default=18051955, help="Random number seed" + ) parser.add_argument( "--noise", type=str, @@ -90,7 +93,10 @@ def cli_parser(): "--ra", type=float, default=0.0, help="Right Ascension (degrees)" ) parser.add_argument( - "--declination", type=float, default=-45.0, help="Declination (degrees)" + "--declination", + type=float, + default=-45.0, + help="Declination (degrees)", ) parser.add_argument( "--nchannels_per_chunk", @@ -118,7 +124,10 @@ def cli_parser(): help="Number of integrations in a chunk to average", ) parser.add_argument( - "--integration_time", type=float, default=0.25, help="Integration time (s)" + "--integration_time", + type=float, + default=0.25, + help="Integration time (s)", ) parser.add_argument( "--time_range", @@ -131,7 +140,8 @@ def cli_parser(): "--input_file", type=str, default="", - help="Full path to the HDF5 file, which contains necessary RFI information for each RFI source.", + help="Full path to the HDF5 file, which contains necessary " + "RFI information for each RFI source.", ) parser.add_argument( "--use_beamgain", @@ -143,31 +153,39 @@ def cli_parser(): "--beamgain_hdf_file", type=str, default="beam_gain_tv_transmitters.hdf5", - help="HDF5 file with beam gain, transmitter, frequency, and pointing (RA, DEC) information.", + help="HDF5 file with beam gain, transmitter, frequency, " + "and pointing (RA, DEC) information.", ) parser.add_argument( "--beamgain_dir", type=str, default="../data/beam_gain", - help="Folder containing multiple Numpy files or the HDF file with beam gain information.", + help="Folder containing multiple Numpy files or the HDF file " + "with beam gain information.", ) parser.add_argument( "--use_antfile", type=str, default="False", - help="Use the antenna file in the rfi data in calculation, otherwise use from RASCIL", + help="Use the antenna file in the rfi data in calculation, " + "otherwise use from RASCIL", ) parser.add_argument( "--antenna_file", type=str, - default="./data/telescope_files/SKA1-LOW_SKO-0000422_Rev3_38m_SKALA4_spot_frequencies.tm/layout_wgs84.txt", + default="./data/telescope_files/" + "SKA1-LOW_SKO-0000422_Rev3_38m_SKALA4_spot_frequencies.tm/" + "layout_wgs84.txt", help="txt file containing antenna locations", ) parser.add_argument( "--write_ms", type=str, default="True", help="Write measurement set?" ) parser.add_argument( - "--msout", type=str, default="./simulate_rfi.ms", help="Name for MeasurementSet" + "--msout", + type=str, + default="./simulate_rfi.ms", + help="Name for MeasurementSet", ) parser.add_argument( "--output_dir", @@ -195,7 +213,8 @@ def add_noise(bvis): :param bvis: BlockVisibility :return: BlockVisibility with noise """ - # The specified sensitivity (effective area / T_sys) is roughly 610 m ^ 2 / K in the range 160 - 200MHz + # The specified sensitivity: + # (effective area / T_sys) is roughly 610 m ^ 2 / K in the range 160 - 200MHz # noqa: E501 # sigma_vis = 2 k T_sys / (area * sqrt(tb)) = 2 k 512 / (610 * sqrt(tb) sens = 610 bt = bvis.channel_bandwidth[0] * bvis.integration_time[0] @@ -223,24 +242,33 @@ def simulate_rfi_image_prop( transmitter_ids=None, rfi_frequencies=None, ): - """Simulate RFI and then average over time and frequency, producing a BlockVisibility + """ + Simulate RFI and then average over time and frequency, + producing a BlockVisibility :param config: RASCIL telescope Configuration e.g. LOW :param times: observation times (hour angles in radians) :param frequency: Frequencies (numpy.array) Hz - :param apparent_transmitter_power: power transmitted by RFI source as seen by an isotropic antenna - (numpy.ndarray) dimensions: [nsource, ntimes, nstations, nfreqs] + :param apparent_transmitter_power: power transmitted by RFI source + as seen by an isotropic antenna + (numpy.ndarray) dimensions: + [nsource, ntimes, nstations, nfreqs] :param channel_bandwidth: Channel bandwidth (numpy.array) Hz :param phasecentre: Phasecentre of observations (as SkyCoord) - :param polarisation_frame: In RASCIL format e.g. PolarisationFrame("stokesI") + :param polarisation_frame: In RASCIL format + e.g. PolarisationFrame("stokesI") :param time_average: Number of integrations to average (int) :param channel_average: Number of channels to average (int) :param noise: Add noise? Add thermal noise (approximate numbers) - :param beam_gain: beam gain (numpy.ndarray or float) to be applied to the apparent power - :param transmitter_coordinates: apparent coordinates of the RFI transmitter (azimuth, elevation, distance) - (numpy.ndarray) dimensions: [nsources, ntimes, nstations, 3] + :param beam_gain: beam gain (numpy.ndarray or float) + to be applied to the apparent power + :param transmitter_coordinates: apparent coordinates of the RFI transmitter + (azimuth, elevation, distance) + (numpy.ndarray) dimensions: + [nsources, ntimes, nstations, 3] :param transmitter_ids: 1D numpy array of the RFI source IDs - :param rfi_frequencies: 1D numpy array of the frequencies where the RFI signal appears + :param rfi_frequencies: 1D numpy array of the frequencies + where the RFI signal appears :return: BlockVisibility containing averaged data """ @@ -249,15 +277,17 @@ def simulate_rfi_image_prop( average_chunks(frequency, numpy.ones_like(frequency), channel_average) )[0] averaged_channel_bandwidth, wts = numpy.array( - average_chunks(channel_bandwidth, numpy.ones_like(frequency), channel_average) + average_chunks( + channel_bandwidth, numpy.ones_like(frequency), channel_average + ) ) averaged_channel_bandwidth *= wts averaged_times = numpy.array( average_chunks(times, numpy.ones_like(times), time_average) )[0] - # times are the hour angles in seconds (!) relative to transit. create_blockvisibility needs - # hour angles in radians + # times are the hour angles in seconds (!) relative to transit. + # create_blockvisibility needs hour angles in radians s2r = numpy.pi / 43200.0 bvis = create_blockvisibility( config, @@ -358,7 +388,8 @@ def load_beam_gain_and_ra_dec(args, source_ids, rfi_freq_chans, log): via user provided args. If beam gain is not used, 1.0 is returned. - Else a numpy array loaded from file, with shape [nrfi_sources x nstations x nchannels]. + Else a numpy array loaded from file, + with shape [nrfi_sources x nstations x nchannels]. If beam gain data are for the station centre, nstations = 1 :param args: user provided arguments @@ -376,7 +407,8 @@ def load_beam_gain_and_ra_dec(args, source_ids, rfi_freq_chans, log): if not args.beamgain_hdf_file == "": hdf_file = f"{bg_direct}{args.beamgain_hdf_file}" log.info( - "Loading beam gain data from HDF file: %s. " "Only loading Stokes I.", + "Loading beam gain data from HDF file: %s. " + "Only loading Stokes I.", hdf_file, ) @@ -391,7 +423,8 @@ def load_beam_gain_and_ra_dec(args, source_ids, rfi_freq_chans, log): if (bg_frequencies != rfi_freq_chans).any(): raise ( ValueError( - "Miss-match between RFI frequencies and beam gain frequencies." + "Miss-match between RFI frequencies " + "and beam gain frequencies." ) ) @@ -406,7 +439,10 @@ def load_beam_gain_and_ra_dec(args, source_ids, rfi_freq_chans, log): beam_gain = numpy.zeros((len(source_ids), 1, len(rfi_freq_chans))) for i, source in enumerate(source_ids): - beam_gain_files = f"{bg_direct}{source}_beam_gain_TIME_SEP_CHAN_SEP_CROSS_POWER_AMP_I_I.txt" + beam_gain_files = ( + f"{bg_direct}{source}_beam_gain_TIME_SEP_CHAN_" + f"SEP_CROSS_POWER_AMP_I_I.txt" + ) beam_gain[i, :, :] = numpy.loadtxt(beam_gain_files) return beam_gain, args.ra, args.declination @@ -417,7 +453,9 @@ def load_beam_gain_and_ra_dec(args, source_ids, rfi_freq_chans, log): return beam_gain, args.ra, args.declination -def get_chunk_start_times(time_range, integration_time, nintegrations_per_chunk, log): +def get_chunk_start_times( + time_range, integration_time, nintegrations_per_chunk, log +): """ Generate start times for each chunk. @@ -433,7 +471,9 @@ def get_chunk_start_times(time_range, integration_time, nintegrations_per_chunk, ) log.info("Start times {}".format(start_times)) - chunk_start_times = [start_times[i : i + 1] for i in range(0, len(start_times))] + chunk_start_times = [ + start_times[i : i + 1] for i in range(0, len(start_times)) + ] log.info("Chunk start times {}".format([c[0] for c in chunk_start_times])) log.info("Processing {nst:d} time chunks".format(nst=len(start_times))) return chunk_start_times @@ -441,7 +481,8 @@ def get_chunk_start_times(time_range, integration_time, nintegrations_per_chunk, def rfi_simulation(args): """ - :param args: user-defined arguments; result of argparse.ArgumentParser.parse_args(). + :param args: user-defined arguments; + result of argparse.ArgumentParser.parse_args() """ settings_file = "old_settings/RFI_propagation_settings.csv" check_for_old_run(vars(args), settings_file) @@ -453,7 +494,9 @@ def rfi_simulation(args): start_epoch = time.asctime() log.info( - "\nSKA LOW RFI simulation using RASCIL\nStarted at {}\n".format(start_epoch) + "\nSKA LOW RFI simulation using RASCIL\nStarted at {}\n".format( + start_epoch + ) ) log.info("Starting LOW low-level RFI simulation") @@ -490,7 +533,8 @@ def rfi_simulation(args): low = generate_ska_antenna_configuration( args.use_antfile, args.antenna_file, rmax, log ) - # xarray update - update all xarray variables to match station skipping in configuration + # xarray update - update all xarray variables + # to match station skipping in configuration sub_low = select_configuration(low, low.names[::station_skip]) phasecentre = SkyCoord( @@ -499,12 +543,13 @@ def rfi_simulation(args): log.info(f"Phasecentre is {phasecentre}") # has to match what came in the HDF5 file - integration_time = args.integration_time # Integration time within a chunk - time_average = ( - args.time_average - ) # Number of integrations to average (??), Integration time after averaging + # Integration time within a chunk + integration_time = args.integration_time + # Number of integrations to average (??), Integration time after averaging + time_average = args.time_average log.info( - "Each chunk has {naint:d} integrations of duration {natime:.2f} (s)".format( + "Each chunk has {naint:d} integrations of " + "duration {natime:.2f} (s)".format( naint=n_time_samples, natime=integration_time ) ) @@ -514,24 +559,32 @@ def rfi_simulation(args): # Frequency and bandwidth related set-up frequency = numpy.linspace( - args.frequency_range[0], args.frequency_range[1], args.nchannels_per_chunk + args.frequency_range[0], + args.frequency_range[1], + args.nchannels_per_chunk, + ) + channel_bandwidth = (frequency[-1] - frequency[0]) / ( + args.nchannels_per_chunk - 1 ) - channel_bandwidth = (frequency[-1] - frequency[0]) / (args.nchannels_per_chunk - 1) channel_average = args.channel_average log.info( - "Each chunk has {narchan:d} frequency channels of width {nabw:.3f} (MHz)".format( + "Each chunk has {narchan:d} frequency channels " + "of width {nabw:.3f} (MHz)".format( narchan=args.nchannels_per_chunk, nabw=channel_bandwidth * 1e-6 ) ) channel_bandwidth = numpy.ones_like(frequency) * channel_bandwidth - # We process the chunks (and accumulate the images) in two stages to avoid a large final reduction + # We process the chunks (and accumulate the images) + # in two stages to avoid a large final reduction vis_graph_list = list() for chunk_start_time in chunk_start_times: for start_time in chunk_start_time: - times = start_time + numpy.arange(n_time_samples) * integration_time + times = ( + start_time + numpy.arange(n_time_samples) * integration_time + ) # Perform the simulation for this chunk vis_graph = rsexecute.execute(simulate_rfi_image_prop)( diff --git a/rfi/rfi_interface/rfi_data_cube.py b/rfi/rfi_interface/rfi_data_cube.py index 153faf63d84f44629aeca0cbcb77b384c7108bc3..08b1ee3018e2a513b4bdd316e9154c5a9832867a 100644 --- a/rfi/rfi_interface/rfi_data_cube.py +++ b/rfi/rfi_interface/rfi_data_cube.py @@ -1,9 +1,9 @@ -import h5py import logging -import numpy as np - from dataclasses import dataclass, field +import h5py +import numpy as np + LOG = logging.getLogger("rfi-logger") @@ -28,7 +28,7 @@ class RFISource: f"{self.__class__.__name__}(" f"id={self.id}, " f"source_type={self.source_type}, " - f"central_frequency={self.central_frequency}[{self._get_unit('central_frequency')}], " + f"central_frequency={self.central_frequency}[{self._get_unit('central_frequency')}], " # noqa: E501 f"bandwidth={self.bandwidth}[{self._get_unit('bandwidth')}], " f"polarisation={self.polarisation}, " f"height={self.height}, " @@ -47,7 +47,8 @@ class RFISignal: metadata={"unit": "Hz"} ) # 1D array of frequency samples - # 3D array: list of data per time per station [ [ [...], [...], [...] ], [ [...], [...], [...] ] ] + # 3D array: list of data per time per station + # [ [ [...], [...], [...] ], [ [...], [...], [...] ] ] apparent_power: np.ndarray = field(metadata={"unit": "W/Hz???"}) # coordinates at SKA stations; 2D-arrays, coordinates per time @@ -55,8 +56,10 @@ class RFISignal: elevation: np.ndarray = field(metadata={"unit": "degree"}) distance: np.ndarray = field(metadata={"unit": "m"}) - # reshape az, el, dist, into a single array: [[[az, el, dist], [az, el, dist]], [[az, el, dist], [az, el, dist]]] - # arr[0][0] --> first time, first station values, arr[0][1] --> first time, second station values, etc. + # reshape az, el, dist, into a single array: + # [[[az, el, dist], [az, el, dist]], [[az, el, dist], [az, el, dist]]] + # arr[0][0] --> first time, first station values, arr[0][1] + # --> first time, second station values, etc. apparent_coordinates: np.ndarray @@ -108,7 +111,12 @@ class DataCube: """ def __init__( - self, times: list, freqs: list, station_ids: list, rmax=None, station_skip=None + self, + times: list, + freqs: list, + station_ids: list, + rmax=None, + station_skip=None, ): self.time_samples = np.array(times).astype("S") self.frequency_channels = np.array(freqs).astype("f") @@ -121,11 +129,15 @@ class DataCube: [None], shape=(1, len(times), len(station_ids), 3), dtype=float ) self.signal = np.empty_like( - [None], shape=(1, len(times), len(station_ids), len(freqs)), dtype=complex + [None], + shape=(1, len(times), len(station_ids), len(freqs)), + dtype=complex, ) self._number_of_sources = 0 - self._rmax = rmax # maximum distance of each station from the array centre + self._rmax = ( + rmax # maximum distance of each station from the array centre + ) self._station_skip = station_skip def validate_input_data(self, input_data): @@ -134,11 +146,15 @@ class DataCube: Data are valid if: - source_id exists - - time samples of the input match the ones that the DataCube was initialized with - - frequency channels of the input match the ones that the DataCube was initialized with - - station ids of the input match the ones that the DataCube was initialized with - - :param input_data: input DataCubePerSource object containing RFI data for a single source + - time samples of the input match the ones that the + DataCube was initialized with + - frequency channels of the input match the ones that + the DataCube was initialized with + - station ids of the input match the ones that + the DataCube was initialized with + + :param input_data: input DataCubePerSource object + containing RFI data for a single source """ if not input_data.rfi_source.id: raise ValueError("New data cube cannot have an empty source_id.") @@ -150,14 +166,18 @@ class DataCube: ) if ( - input_data.rfi_signal.frequency.astype("f") != self.frequency_channels + input_data.rfi_signal.frequency.astype("f") + != self.frequency_channels ).all(): raise ValueError( - "The frequency_channels that initialized the DataCube object are " - "different from the appended DataCubePerSource() frequencies." + "The frequency_channels that initialized " + "the DataCube object are different from the " + "appended DataCubePerSource() frequencies." ) - if (input_data.rfi_signal.station_id.astype("S") != self.station_id).all(): + if ( + input_data.rfi_signal.station_id.astype("S") != self.station_id + ).all(): raise ValueError( "The station ids that initialized the DataCube object are " "different from the appended DataCubePerSource() station ids." @@ -167,7 +187,8 @@ class DataCube: """ Append data from a DataCubePerSource object to the existing arrays. - :param new_rfi_data: input DataCubePerSource object containing RFI data for a single source + :param new_rfi_data: input DataCubePerSource object + containing RFI data for a single source """ self.validate_input_data(new_rfi_data) @@ -190,8 +211,11 @@ class DataCube: ] = new_rfi_data.rfi_signal.apparent_coordinates.astype("f") except IndexError: - # when a new source is added, the index of the array has to be increased by one first - self.coordinates = np.append(self.coordinates, self.coordinates, axis=0) + # when a new source is added, the index of the + # array has to be increased by one first + self.coordinates = np.append( + self.coordinates, self.coordinates, axis=0 + ) self.signal = np.append(self.signal, self.signal, axis=0) self.signal[ @@ -219,7 +243,9 @@ class DataCube: hdf5_file["signal"].attrs["dimensions"] = [ "Nsources x Ntimes x Nstations x Nfreqs" ] - hdf5_file["station_id"].attrs["max_station_distance_from_centre"] = self._rmax + hdf5_file["station_id"].attrs[ + "max_station_distance_from_centre" + ] = self._rmax hdf5_file["station_id"].attrs["station_skip"] = self._station_skip def export_to_hdf5(self, filename): @@ -281,7 +307,8 @@ class BeamGainDataCube: if self._beam_gain is not None: LOG.warning( - f"Overwriting existing data in {self.__class__.__name__}.beam_gain." + f"Overwriting existing data in " + f"{self.__class__.__name__}.beam_gain." ) cont = _get_input() if cont.lower() != "yes": @@ -297,8 +324,8 @@ class BeamGainDataCube: ): raise ValueError( "Input data has the wrong dimensions. It needs to match: \n" - f"{len(self.rfi_source_ids), len(self.freq_channels), self.num_stations, 4} \n" - f"[Nrfi_sources x Nska_stations x Nfreq_channels x 4 (Nstokes)]" + f"{len(self.rfi_source_ids), len(self.freq_channels), self.num_stations, 4} \n" # noqa: E501 + f"[Nrfi_sources x Nska_stations x Nfreq_channels x 4 (Nstokes)]" # noqa: E501 ) self._beam_gain = value @@ -308,7 +335,8 @@ class BeamGainDataCube: hdf5_file[k.lstrip("_")] = self.__dict__[k] hdf5_file["beam_gain"].attrs["dimensions"] = [ - "Nrfi_sources x Nska_stations x Nfreq_channels x 4 (stokesI, stokesQ, stokesU, stokesV)" + "Nrfi_sources x Nska_stations x Nfreq_channels " + "x 4 (stokesI, stokesQ, stokesU, stokesV)" ] def export_to_hdf5(self, filename): diff --git a/rfi/rfi_interface/rfi_source_signal_interface.py b/rfi/rfi_interface/rfi_source_signal_interface.py index 927277e163df2ce899961069514b7ea1f809b052..0cede063ae3c858d0197f30d242eec36f1c9a1b1 100644 --- a/rfi/rfi_interface/rfi_source_signal_interface.py +++ b/rfi/rfi_interface/rfi_source_signal_interface.py @@ -20,12 +20,16 @@ Options: <freq_start> Start of Frequency range [MHz]. Default: 170.5 <freq_end> End of Frequency range [MHz]. Default: 184.5 -""" +""" # noqa: E501 import sys + import numpy as np from docopt import docopt -from rfi.pycraf_scripts.SKA_low_RFI_propagation import rfi_attenuation, cli_parser +from rfi.pycraf_scripts.SKA_low_RFI_propagation import ( + cli_parser, + rfi_attenuation, +) from rfi.rfi_interface.rfi_data_cube import DataCube N_CHANNELS = "n_channels" @@ -58,7 +62,11 @@ def rfi_interface(source_type, **kwargs): if source_type == "tv_antenna": args_list = [] - [args_list.extend([k, v]) for k, v in kwargs.items() if k != FREQ_RANGE] + [ + args_list.extend([k, v]) + for k, v in kwargs.items() + if k != FREQ_RANGE + ] args_list.extend( [ FREQ_RANGE, @@ -78,7 +86,11 @@ def rfi_interface(source_type, **kwargs): station_ids = apparent_power_results[0].rfi_signal.station_id times = [None] * int(args.n_time_chunks) cube = DataCube( - times, freqs, station_ids, rmax=args.rmax, station_skip=args.station_skip + times, + freqs, + station_ids, + rmax=args.rmax, + station_skip=args.station_skip, ) for source_data in apparent_power_results: @@ -86,7 +98,9 @@ def rfi_interface(source_type, **kwargs): cube.append_data(source_data) elif source_type == "aircraft": - raise NotImplementedError("Aircraft RFI calculations are not implemented yet.") + raise NotImplementedError( + "Aircraft RFI calculations are not implemented yet." + ) return cube @@ -104,7 +118,9 @@ def export_data(attenuation_data, source_type): def validate_frequency(arguments): validate = [arguments[k] for k in FREQ_KEYS if arguments[k] is not None] if len(validate) > 0 and len(validate) != 3: - raise ValueError(f"If you provide one of {FREQ_KEYS}, you need to provide all.") + raise ValueError( + f"If you provide one of {FREQ_KEYS}, you need to provide all." + ) if not all([arguments[k] for k in FREQ_KEYS]): # default values for SKA Low in SKA_low_RFI_propagation.py @@ -113,7 +129,10 @@ def validate_frequency(arguments): else: arguments[f"--{N_CHANNELS}"] = str(arguments[f"<{N_CHANNELS}>"]) - arguments[FREQ_RANGE] = [str(arguments[FREQ_START]), str(arguments[FREQ_END])] + arguments[FREQ_RANGE] = [ + str(arguments[FREQ_START]), + str(arguments[FREQ_END]), + ] # remove keys that cannot be used by RFI code arguments.pop(f"<{N_CHANNELS}>") @@ -124,9 +143,12 @@ def validate_frequency(arguments): def main(argv): args = docopt(__doc__, argv=argv[1:], help=True) - source_type = [x for x in ACCEPTED_SOURCES if args[x]][0] # only one can be true + source_type = [x for x in ACCEPTED_SOURCES if args[x]][ + 0 + ] # only one can be true - # remove keys that will not be an input to any of the rfi_attenuation functions + # remove keys that will not be an input to any + # of the rfi_attenuation functions args.pop("--help", None) args.pop("-h", None) for source in ACCEPTED_SOURCES: diff --git a/scripts/animate_ms.py b/scripts/animate_ms.py index e2b3e86b96f6c0249f9002c3ea0ffbb334eea0aa..040f131c82b9932f5c00823d76684f8ad63c2c7a 100755 --- a/scripts/animate_ms.py +++ b/scripts/animate_ms.py @@ -30,7 +30,9 @@ class Plotter: self._num_frames = 0 self._title = "" - def animate(self, imager_settings, ms_names, title="", fps=10, filename="out.mp4"): + def animate( + self, imager_settings, ms_names, title="", fps=10, filename="out.mp4" + ): """Function to generate the animation. Args: @@ -144,7 +146,9 @@ class Plotter: if start_row >= num_rows or start_row + num_baselines > num_rows: continue (u, v, w) = self._ms_list[i].read_coords(start_row, num_baselines) - vis = self._ms_list[i].read_column("DATA", start_row, num_baselines) + vis = self._ms_list[i].read_column( + "DATA", start_row, num_baselines + ) num_pols = vis.shape[-1] # Create settings for the imager. @@ -159,7 +163,9 @@ class Plotter: ) imager = oskar.Imager(settings=settings) imager.set_vis_frequency(freq_start_hz, freq_inc_hz, num_channels) - imager.update(u, v, w, vis, end_channel=num_channels - 1, num_pols=num_pols) + imager.update( + u, v, w, vis, end_channel=num_channels - 1, num_pols=num_pols + ) data = imager.finalise(return_images=1) # Update the plot panel and colourbar. @@ -177,7 +183,10 @@ def main(): "ms_names", metavar="MS", nargs="+", help="Measurement Set path(s)" ) parser.add_argument( - "--fov_deg", type=float, default=0.5, help="Field of view to image, in degrees" + "--fov_deg", + type=float, + default=0.5, + help="Field of view to image, in degrees", ) parser.add_argument( "--size", type=int, default=256, help="Image side length, in pixels" @@ -194,7 +203,9 @@ def main(): # Make animation. plotter = Plotter() - plotter.animate(imager_settings, args.ms_names, args.title, args.fps, args.out) + plotter.animate( + imager_settings, args.ms_names, args.title, args.fps, args.out + ) if __name__ == "__main__": diff --git a/scripts/create_ionospheric_screen.py b/scripts/create_ionospheric_screen.py index c010a03f55212901d993154ce81a0a1b20fbab1d..befd25bb843ff871e485360bc602d0eadcce8c30 100644 --- a/scripts/create_ionospheric_screen.py +++ b/scripts/create_ionospheric_screen.py @@ -14,13 +14,22 @@ def main(): # Get command line arguments. parser = argparse.ArgumentParser() parser.add_argument( - "--num_times", type=int, default=240, help="number of time samples to generate" + "--num_times", + type=int, + default=240, + help="number of time samples to generate", ) parser.add_argument( - "--interval_sec", type=float, default=60.0, help="interval between time samples" + "--interval_sec", + type=float, + default=60.0, + help="interval between time samples", ) parser.add_argument( - "--pixel_size_m", type=float, default=100.0, help="pixel size, in metres" + "--pixel_size_m", + type=float, + default=100.0, + help="pixel size, in metres", ) parser.add_argument( "--screen_width_km", @@ -29,7 +38,10 @@ def main(): help="width of screen, in kilometres", ) parser.add_argument( - "--scale_size_km", type=float, default=5.0, help="scale size, in kilometres" + "--scale_size_km", + type=float, + default=5.0, + help="scale size, in kilometres", ) parser.add_argument("filename", help="name of FITS file to write") args = parser.parse_args() @@ -39,7 +51,9 @@ def main(): bmax = 0.1 * screen_width_metres # sub-aperture(?) size. sampling = args.pixel_size_m m = int(bmax / sampling) # Pixels per sub-aperture (200). - n = int(screen_width_metres / bmax) # Sub-apertures across the screen (10). + n = int( + screen_width_metres / bmax + ) # Sub-apertures across the screen (10). num_pix = n * m pscale = screen_width_metres / (n * m) # Pixel scale (100 m/pixel). print("Number of pixels %d, pixel size %.3f m" % (num_pix, pscale)) diff --git a/scripts/generate_station_layouts.py b/scripts/generate_station_layouts.py index 4cdb6be9f31062a94d43934ffc682bac37d2d86a..be8d8dbbfb55f825ea38b2784e6387a9c0ac6f89 100644 --- a/scripts/generate_station_layouts.py +++ b/scripts/generate_station_layouts.py @@ -43,7 +43,7 @@ def random_layout(n, r_max, min_sep, timeout=None): # Get trial coordinates. xt, yt = tuple(numpy.random.rand(2) * 2.0 * r_max - r_max) - if ((xt ** 2 + yt ** 2) ** 0.5 + min_sep / 2.0) > r_max: + if ((xt**2 + yt**2) ** 0.5 + min_sep / 2.0) > r_max: continue # Get tile indices and tile ranges to check. @@ -63,7 +63,7 @@ def random_layout(n, r_max, min_sep, timeout=None): for _ in range(grid_count[ky, kx]): dx = xt - xy[i_other, 0] dy = yt - xy[i_other, 1] - d_other = (dx ** 2 + dy ** 2) ** 0.5 + d_other = (dx**2 + dy**2) ** 0.5 d_min = min(d_min, d_other) i_other = grid_next[i_other] diff --git a/scripts/low_dd_effects_validation.py b/scripts/low_dd_effects_validation.py index 468838ec7a3395e7414e4d2d4ab1d552085f4c86..71dd64d056015b9fdda52a590a65829f27ab4b74 100644 --- a/scripts/low_dd_effects_validation.py +++ b/scripts/low_dd_effects_validation.py @@ -63,7 +63,9 @@ class Plotter: self._ms1 = ms1 def animate(self, num_times, fps, filename): - fig, axes = plt.subplots(nrows=len(self._freqs), ncols=4, figsize=(16, 8)) + fig, axes = plt.subplots( + nrows=len(self._freqs), ncols=4, figsize=(16, 8) + ) self._fig = fig self._axes = axes.flatten() anim = animation.FuncAnimation( @@ -157,15 +159,27 @@ class Plotter: # Make the image. print("Generating frame %d, panel %d" % (frame, panel)) imager = oskar.Imager(settings=settings) - imager.set_vis_frequency(freq_start_hz, freq_inc_hz, num_channels) + imager.set_vis_frequency( + freq_start_hz, freq_inc_hz, num_channels + ) imager.set_vis_phase_centre(ra0_deg, dec0_deg) imager.coords_only = True imager.update( - u, v, w, vis[i], end_channel=num_channels - 1, num_pols=4 + u, + v, + w, + vis[i], + end_channel=num_channels - 1, + num_pols=4, ) imager.coords_only = False imager.update( - u, v, w, vis[i], end_channel=num_channels - 1, num_pols=4 + u, + v, + w, + vis[i], + end_channel=num_channels - 1, + num_pols=4, ) data = imager.finalise(return_images=1) @@ -222,8 +236,14 @@ def main(): for obs in observations: for field in fields: # Open the Measurement Sets. - ms0_name = "SKA_LOW_SIM_%s_%s_ionosphere_off_GLEAM.MS" % (obs, field) - ms1_name = "SKA_LOW_SIM_%s_%s_ionosphere_on_GLEAM.MS" % (obs, field) + ms0_name = "SKA_LOW_SIM_%s_%s_ionosphere_off_GLEAM.MS" % ( + obs, + field, + ) + ms1_name = "SKA_LOW_SIM_%s_%s_ionosphere_on_GLEAM.MS" % ( + obs, + field, + ) ms0 = oskar.MeasurementSet.open(ms0_name, readonly=True) ms1 = oskar.MeasurementSet.open(ms1_name, readonly=True) @@ -237,7 +257,9 @@ def main(): # Make animation. weighting = base_settings["image/weighting"] filename = "SKA_LOW_SIM_%s_%s_%s.mp4" % (obs, field, weighting) - plotter = Plotter(base_settings, obs, field, freqs, sources, ms0, ms1) + plotter = Plotter( + base_settings, obs, field, freqs, sources, ms0, ms1 + ) plotter.animate(num_times, 10, filename) diff --git a/scripts/plot_beams.py b/scripts/plot_beams.py index 93e8435810eee2067738ad1e6327a3728c6407cb..cc7bbfa2f8b3834d76ee468618989aac2aeaf6c3 100644 --- a/scripts/plot_beams.py +++ b/scripts/plot_beams.py @@ -25,7 +25,9 @@ def plot_panel(ax, image, title, cmap): """Plots a single panel.""" im = ax.imshow(numpy.squeeze(image), cmap=cmap) plt.colorbar(im, format="%.2e") - plt.tick_params(labelcolor="none", top="off", bottom="off", left="off", right="off") + plt.tick_params( + labelcolor="none", top="off", bottom="off", left="off", right="off" + ) ax.set_xticks([]) ax.set_yticks([]) ax.invert_yaxis() @@ -76,7 +78,9 @@ def make_plots(title, glob_pattern, out_basename): plt.close("all") -def run_single(current_settings, tel_name, freq, pointing_name, pointing_config): +def run_single( + current_settings, tel_name, freq, pointing_name, pointing_config +): """Runs for a single configuration.""" # Name of the application to run, and a settings file for it. app = "oskar_sim_beam_pattern" @@ -110,7 +114,9 @@ def run_single(current_settings, tel_name, freq, pointing_name, pointing_config) if pointing_name: title += " (%s)" % pointing_name make_plots( - title=title, glob_pattern=root_path + "*_AMP*", out_basename=root_path + "_amp" + title=title, + glob_pattern=root_path + "*_AMP*", + out_basename=root_path + "_amp", ) make_plots( title=title, @@ -220,7 +226,11 @@ def main(): # Station beam: loop over pointing directions. for pointing_name, pointing_config in pointings.items(): run_single( - current_settings, tel_name, freq, pointing_name, pointing_config + current_settings, + tel_name, + freq, + pointing_name, + pointing_config, ) diff --git a/scripts/plot_stokes.py b/scripts/plot_stokes.py index ef2bcc55ab741072fc64697f49be4e7a04f823b3..4b1eb5196ea950a843978ad23aae4c7ed032e0d9 100644 --- a/scripts/plot_stokes.py +++ b/scripts/plot_stokes.py @@ -13,7 +13,9 @@ def plot_panel(ax, image, title, cmap): """Plots a single panel.""" im = ax.imshow(numpy.squeeze(image), cmap=cmap) plt.colorbar(im, format="%.2e") - plt.tick_params(labelcolor="none", top="off", bottom="off", left="off", right="off") + plt.tick_params( + labelcolor="none", top="off", bottom="off", left="off", right="off" + ) ax.set_xticks([]) ax.set_yticks([]) ax.invert_yaxis() @@ -55,7 +57,9 @@ def main(): # Common axis labels. ax = fig.add_subplot(111, frameon=False) - plt.tick_params(labelcolor="none", top="off", bottom="off", left="off", right="off") + plt.tick_params( + labelcolor="none", top="off", bottom="off", left="off", right="off" + ) # Hide ticks for this new axis. ax.set_xticks([]) ax.set_yticks([]) diff --git a/scripts/power_spectrum_test.py b/scripts/power_spectrum_test.py index feaaca9cf0abbe7aab3c04d90dd1b0bd6ead697d..291f32c26f978764c2ae7538d2b818323c3243e4 100644 --- a/scripts/power_spectrum_test.py +++ b/scripts/power_spectrum_test.py @@ -20,7 +20,7 @@ def main(file_name): image = fits_file[0].data[0] cellsize = abs(fits_file[0].header["CDELT1"]) * numpy.pi / 180.0 spectrum = numpy.abs(numpy.fft.fftshift(numpy.fft.fft2(image))) - profile = radial_profile(spectrum ** 2) + profile = radial_profile(spectrum**2) # Delta_u = 1 / (N * Delta_theta) cellsize_uv = 1.0 / (image.shape[0] * cellsize) lambda_max = cellsize_uv * len(profile) @@ -40,5 +40,7 @@ def main(file_name): if __name__ == "__main__": if len(sys.argv) < 2: - raise RuntimeError("Need path of input FITS image file on command line.") + raise RuntimeError( + "Need path of input FITS image file on command line." + ) main(sys.argv[1]) diff --git a/scripts/residual_image_simulator.py b/scripts/residual_image_simulator.py index cdf174aa4538a9ba70b685d16fdfe1b71b26725e..523fd79b5f78419e30661855eb05507ade70ebf5 100644 --- a/scripts/residual_image_simulator.py +++ b/scripts/residual_image_simulator.py @@ -144,10 +144,14 @@ class ResidualImageSimulator(oskar.Interferometer): if not self.coords_only: if self._ref_vis: # Read the reference visibility data block. - self._ref_block.read(self._ref_hdr, self._ref_handle, block_index) + self._ref_block.read( + self._ref_hdr, self._ref_handle, block_index + ) # Subtract reference block from this one. - block.cross_correlations()[:] -= self._ref_block.cross_correlations()[:] + block.cross_correlations()[ + : + ] -= self._ref_block.cross_correlations()[:] # Write (residual) block to any open files. self.write_block(block, block_index) diff --git a/scripts/run_sims_beam_errors.py b/scripts/run_sims_beam_errors.py index 2a11d2a020808f58cafa4e4d63dcd0763f818f17..2039b2eb093dd761686ff2e590c09d9e4d032bf3 100644 --- a/scripts/run_sims_beam_errors.py +++ b/scripts/run_sims_beam_errors.py @@ -30,15 +30,132 @@ def append_bright_sources(sky): numpy.array( ( [266.41683, -29.00781, 2000, 0, 0, 0, 0, 0, 0, 3600, 3600, 0], - [50.67375, -37.20833, 528, 0, 0, 0, 178e6, -0.51, 0, 0, 0, 0], # For - [201.36667, -43.01917, 1370, 0, 0, 0, 200e6, -0.50, 0, 0, 0, 0], # Cen - [139.52500, -12.09556, 280, 0, 0, 0, 200e6, -0.96, 0, 0, 0, 0], # Hyd - [79.95833, -45.77889, 390, 0, 0, 0, 200e6, -0.99, 0, 0, 0, 0], # Pic - [252.78333, 4.99250, 377, 0, 0, 0, 200e6, -1.07, 0, 0, 0, 0], # Her - [187.70417, 12.39111, 861, 0, 0, 0, 200e6, -0.86, 0, 0, 0, 0], # Vir - [83.63333, 22.01444, 1340, 0, 0, 0, 200e6, -0.22, 0, 0, 0, 0], # Tau - [299.86667, 40.73389, 7920, 0, 0, 0, 200e6, -0.78, 0, 0, 0, 0], # Cyg - [350.86667, 58.81167, 11900, 0, 0, 0, 200e6, -0.41, 0, 0, 0, 0], # Cas + [ + 50.67375, + -37.20833, + 528, + 0, + 0, + 0, + 178e6, + -0.51, + 0, + 0, + 0, + 0, + ], # For + [ + 201.36667, + -43.01917, + 1370, + 0, + 0, + 0, + 200e6, + -0.50, + 0, + 0, + 0, + 0, + ], # Cen + [ + 139.52500, + -12.09556, + 280, + 0, + 0, + 0, + 200e6, + -0.96, + 0, + 0, + 0, + 0, + ], # Hyd + [ + 79.95833, + -45.77889, + 390, + 0, + 0, + 0, + 200e6, + -0.99, + 0, + 0, + 0, + 0, + ], # Pic + [ + 252.78333, + 4.99250, + 377, + 0, + 0, + 0, + 200e6, + -1.07, + 0, + 0, + 0, + 0, + ], # Her + [ + 187.70417, + 12.39111, + 861, + 0, + 0, + 0, + 200e6, + -0.86, + 0, + 0, + 0, + 0, + ], # Vir + [ + 83.63333, + 22.01444, + 1340, + 0, + 0, + 0, + 200e6, + -0.22, + 0, + 0, + 0, + 0, + ], # Tau + [ + 299.86667, + 40.73389, + 7920, + 0, + 0, + 0, + 200e6, + -0.78, + 0, + 0, + 0, + 0, + ], # Cyg + [ + 350.86667, + 58.81167, + 11900, + 0, + 0, + 0, + 200e6, + -0.41, + 0, + 0, + 0, + 0, + ], # Cas ) ) ) @@ -61,7 +178,9 @@ def make_vis_data(settings, sky, tel): sim.run() -def make_diff_image_stats(filename1, filename2, use_w_projection, out_image_root=None): +def make_diff_image_stats( + filename1, filename2, use_w_projection, out_image_root=None +): """Make an image of the difference between two visibility data sets. This function assumes that the observation parameters for both data sets @@ -113,10 +232,10 @@ def make_diff_image_stats(filename1, filename2, use_w_projection, out_image_root "image_medianabs": numpy.median(numpy.abs(image)), "image_mean": numpy.mean(image), "image_std": numpy.std(image), - "image_rms": numpy.sqrt(numpy.mean(image ** 2)), + "image_rms": numpy.sqrt(numpy.mean(image**2)), "image_centre_mean": numpy.mean(centre), "image_centre_std": numpy.std(centre), - "image_centre_rms": numpy.sqrt(numpy.mean(centre ** 2)), + "image_centre_rms": numpy.sqrt(numpy.mean(centre**2)), } @@ -169,7 +288,12 @@ def make_plot( # Scatter overlay. sp = ax1.scatter( - special_gain, special_phase, c=special_data, cmap="plasma", norm=norm, zorder=10 + special_gain, + special_phase, + c=special_data, + cmap="plasma", + norm=norm, + zorder=10, ) for i, val in enumerate(special_data): ax1.annotate( @@ -208,7 +332,14 @@ def make_plot( def run_single( - prefix_field, settings, sky, tel, gain_std_dB, phase_std_deg, out0_name, results + prefix_field, + settings, + sky, + tel, + gain_std_dB, + phase_std_deg, + out0_name, + results, ): """Run a single simulation and generate image statistics for it.""" out = "%s_%.3f_dB_%.2f_deg" % (prefix_field, gain_std_dB, phase_std_deg) @@ -268,7 +399,9 @@ def run_set( ra_deg = float(settings["observation/phase_centre_ra_deg"]) dec_deg = float(settings["observation/phase_centre_dec_deg"]) length_sec = float(settings["observation/length"]) - settings["observation/start_time_utc"] = get_start_time(ra_deg, length_sec) + settings["observation/start_time_utc"] = get_start_time( + ra_deg, length_sec + ) tel.set_phase_centre(ra_deg, dec_deg) # Load or create the sky model. diff --git a/scripts/run_sims_element_patterns.py b/scripts/run_sims_element_patterns.py index e2a807c1ee25d4453c6065ed13f74f55309b7664..a6c31dab3db4bbcbbb3b6d28cee0c531772f4200 100644 --- a/scripts/run_sims_element_patterns.py +++ b/scripts/run_sims_element_patterns.py @@ -31,15 +31,132 @@ def bright_sources(): return numpy.array( ( [266.41683, -29.00781, 2000, 0, 0, 0, 0, 0, 0, 3600, 3600, 0], - [50.67375, -37.20833, 528, 0, 0, 0, 178e6, -0.51, 0, 0, 0, 0], # For - [201.36667, -43.01917, 1370, 0, 0, 0, 200e6, -0.50, 0, 0, 0, 0], # Cen - [139.52500, -12.09556, 280, 0, 0, 0, 200e6, -0.96, 0, 0, 0, 0], # Hyd - [79.95833, -45.77889, 390, 0, 0, 0, 200e6, -0.99, 0, 0, 0, 0], # Pic - [252.78333, 4.99250, 377, 0, 0, 0, 200e6, -1.07, 0, 0, 0, 0], # Her - [187.70417, 12.39111, 861, 0, 0, 0, 200e6, -0.86, 0, 0, 0, 0], # Vir - [83.63333, 22.01444, 1340, 0, 0, 0, 200e6, -0.22, 0, 0, 0, 0], # Tau - [299.86667, 40.73389, 7920, 0, 0, 0, 200e6, -0.78, 0, 0, 0, 0], # Cyg - [350.86667, 58.81167, 11900, 0, 0, 0, 200e6, -0.41, 0, 0, 0, 0], # Cas + [ + 50.67375, + -37.20833, + 528, + 0, + 0, + 0, + 178e6, + -0.51, + 0, + 0, + 0, + 0, + ], # For + [ + 201.36667, + -43.01917, + 1370, + 0, + 0, + 0, + 200e6, + -0.50, + 0, + 0, + 0, + 0, + ], # Cen + [ + 139.52500, + -12.09556, + 280, + 0, + 0, + 0, + 200e6, + -0.96, + 0, + 0, + 0, + 0, + ], # Hyd + [ + 79.95833, + -45.77889, + 390, + 0, + 0, + 0, + 200e6, + -0.99, + 0, + 0, + 0, + 0, + ], # Pic + [ + 252.78333, + 4.99250, + 377, + 0, + 0, + 0, + 200e6, + -1.07, + 0, + 0, + 0, + 0, + ], # Her + [ + 187.70417, + 12.39111, + 861, + 0, + 0, + 0, + 200e6, + -0.86, + 0, + 0, + 0, + 0, + ], # Vir + [ + 83.63333, + 22.01444, + 1340, + 0, + 0, + 0, + 200e6, + -0.22, + 0, + 0, + 0, + 0, + ], # Tau + [ + 299.86667, + 40.73389, + 7920, + 0, + 0, + 0, + 200e6, + -0.78, + 0, + 0, + 0, + 0, + ], # Cyg + [ + 350.86667, + 58.81167, + 11900, + 0, + 0, + 0, + 200e6, + -0.41, + 0, + 0, + 0, + 0, + ], # Cas ) ) @@ -81,11 +198,15 @@ def make_sky_model(sky0, settings, radius_deg, flux_min_outer_jy): sky_outer.num_sources, ) sky_outer.append(sky_inner) - LOG.info("Number of sources in output sky model: %d", sky_outer.num_sources) + LOG.info( + "Number of sources in output sky model: %d", sky_outer.num_sources + ) return sky_outer -def make_plot(prefix, field_name, metric_key, results, axis_length, axis_bandwidth): +def make_plot( + prefix, field_name, metric_key, results, axis_length, axis_bandwidth +): """Plot selected results.""" # Axis setup. ax1 = plt.subplot(111) @@ -143,13 +264,19 @@ def create_imagers_multi_channel( channel_end = channel_start + channel_step - 1 if channel_end >= num_channels - 1: channel_end = num_channels - 1 - key = "length_%.0f_channels_%d-%d" % (length_sec, channel_start, channel_end) + key = "length_%.0f_channels_%d-%d" % ( + length_sec, + channel_start, + channel_end, + ) if key in results: continue imager = oskar.Imager(settings=settings_img) freq_min = freq_start + freq_inc * (channel_start - 0.5) freq_max = freq_start + freq_inc * (channel_end + 0.5) - LOG.info('"%s" has frequency range %.5e to %.5e Hz', key, freq_min, freq_max) + LOG.info( + '"%s" has frequency range %.5e to %.5e Hz', key, freq_min, freq_max + ) imager.channel_snapshots = True imager.grid_on_gpu = True imager.freq_min_hz = freq_min @@ -158,7 +285,9 @@ def create_imagers_multi_channel( imager.time_max_utc = time_max imager.output_root = prefix_field + "_" + key imagers.append((key, imager)) - LOG.info("Will make %d images of %d channels each.", len(imagers), channel_step) + LOG.info( + "Will make %d images of %d channels each.", len(imagers), channel_step + ) return imagers @@ -198,7 +327,9 @@ def read_coords(out0_name, imagers): # Read both visibility blocks. LOG.info("Reading blocks %d/%d", i_block + 1, hdr1.num_blocks) tasks_read.append( - executor.submit(block_diff[i_block % 2].read, hdr1, handle1, i_block) + executor.submit( + block_diff[i_block % 2].read, hdr1, handle1, i_block + ) ) # Ensure reads have finished. @@ -269,7 +400,13 @@ def read_data(out0_name, out_name, prefix_field, imagers): def run_set_new( - base_settings_sim, base_settings_img, sky0, prefix_field, field, length_sec, results + base_settings_sim, + base_settings_img, + sky0, + prefix_field, + field, + length_sec, + results, ): """Runs a set of simulations for one field.""" # Create settings. @@ -415,7 +552,10 @@ def run_set( freq_min = freq_mid - (bandwidth * 1e6 + freq_inc) * 0.5 freq_max = freq_mid + (bandwidth * 1e6 + freq_inc) * 0.5 LOG.info( - '"%s" has frequency range %.5e to %.5e Hz', key, freq_min, freq_max + '"%s" has frequency range %.5e to %.5e Hz', + key, + freq_min, + freq_max, ) imager.freq_min_hz = freq_min imager.freq_max_hz = freq_max @@ -442,8 +582,12 @@ def run_set( block1.read(hdr1, handle1, i_block) tasks = [] for i, (key, imager) in enumerate(imagers): - LOG.info('Updating imager %d/%d "%s"...', i + 1, len(imagers), key) - tasks.append(executor.submit(imager.update_from_block, hdr1, block1)) + LOG.info( + 'Updating imager %d/%d "%s"...', i + 1, len(imagers), key + ) + tasks.append( + executor.submit(imager.update_from_block, hdr1, block1) + ) concurrent.futures.wait(tasks) for i, (key, imager) in enumerate(imagers): imager.coords_only = False @@ -465,7 +609,9 @@ def run_set( block1.cross_correlations()[...] -= block2.cross_correlations() tasks = [] for _, imager in imagers: - tasks.append(executor.submit(imager.update_from_block, hdr1, block1)) + tasks.append( + executor.submit(imager.update_from_block, hdr1, block1) + ) concurrent.futures.wait(tasks) executor.shutdown() del handle1, handle2, hdr1, hdr2, block1, block2 @@ -487,10 +633,10 @@ def run_set( "image_medianabs": numpy.median(numpy.abs(image)), "image_mean": numpy.mean(image), "image_std": numpy.std(image), - "image_rms": numpy.sqrt(numpy.mean(image ** 2)), + "image_rms": numpy.sqrt(numpy.mean(image**2)), "image_centre_mean": numpy.mean(centre), "image_centre_std": numpy.std(centre), - "image_centre_rms": numpy.sqrt(numpy.mean(centre ** 2)), + "image_centre_rms": numpy.sqrt(numpy.mean(centre**2)), } del image, output @@ -528,7 +674,9 @@ def main(): } base_settings_img = { "image": { - "double_precision": base_settings_sim["simulator"]["double_precision"], + "double_precision": base_settings_sim["simulator"][ + "double_precision" + ], "size": 8192, "fov_deg": 5.0, "fft/use_gpu": True, @@ -581,7 +729,8 @@ def main(): prefix_field = prefix + "_" + field_name results = {} json_file = ( - prefix_field + "_results_no_errors_SKALA4_coupled_SKALA4_single.json" + prefix_field + + "_results_no_errors_SKALA4_coupled_SKALA4_single.json" ) if os.path.exists(json_file): with open(json_file, "r") as input_file: diff --git a/scripts/run_sims_ionosphere_spot_freqs.py b/scripts/run_sims_ionosphere_spot_freqs.py index f05b3d6e2e924e21810e4b3ce49cdaf4788af7c4..5f98b5f43faeccc6576ba3e84bde3cbd9109249b 100644 --- a/scripts/run_sims_ionosphere_spot_freqs.py +++ b/scripts/run_sims_ionosphere_spot_freqs.py @@ -31,15 +31,132 @@ def bright_sources(): return numpy.array( ( [266.41683, -29.00781, 2000, 0, 0, 0, 0, 0, 0, 3600, 3600, 0], - [50.67375, -37.20833, 528, 0, 0, 0, 178e6, -0.51, 0, 0, 0, 0], # For - [201.36667, -43.01917, 1370, 0, 0, 0, 200e6, -0.50, 0, 0, 0, 0], # Cen - [139.52500, -12.09556, 280, 0, 0, 0, 200e6, -0.96, 0, 0, 0, 0], # Hyd - [79.95833, -45.77889, 390, 0, 0, 0, 200e6, -0.99, 0, 0, 0, 0], # Pic - [252.78333, 4.99250, 377, 0, 0, 0, 200e6, -1.07, 0, 0, 0, 0], # Her - [187.70417, 12.39111, 861, 0, 0, 0, 200e6, -0.86, 0, 0, 0, 0], # Vir - [83.63333, 22.01444, 1340, 0, 0, 0, 200e6, -0.22, 0, 0, 0, 0], # Tau - [299.86667, 40.73389, 7920, 0, 0, 0, 200e6, -0.78, 0, 0, 0, 0], # Cyg - [350.86667, 58.81167, 11900, 0, 0, 0, 200e6, -0.41, 0, 0, 0, 0], # Cas + [ + 50.67375, + -37.20833, + 528, + 0, + 0, + 0, + 178e6, + -0.51, + 0, + 0, + 0, + 0, + ], # For + [ + 201.36667, + -43.01917, + 1370, + 0, + 0, + 0, + 200e6, + -0.50, + 0, + 0, + 0, + 0, + ], # Cen + [ + 139.52500, + -12.09556, + 280, + 0, + 0, + 0, + 200e6, + -0.96, + 0, + 0, + 0, + 0, + ], # Hyd + [ + 79.95833, + -45.77889, + 390, + 0, + 0, + 0, + 200e6, + -0.99, + 0, + 0, + 0, + 0, + ], # Pic + [ + 252.78333, + 4.99250, + 377, + 0, + 0, + 0, + 200e6, + -1.07, + 0, + 0, + 0, + 0, + ], # Her + [ + 187.70417, + 12.39111, + 861, + 0, + 0, + 0, + 200e6, + -0.86, + 0, + 0, + 0, + 0, + ], # Vir + [ + 83.63333, + 22.01444, + 1340, + 0, + 0, + 0, + 200e6, + -0.22, + 0, + 0, + 0, + 0, + ], # Tau + [ + 299.86667, + 40.73389, + 7920, + 0, + 0, + 0, + 200e6, + -0.78, + 0, + 0, + 0, + 0, + ], # Cyg + [ + 350.86667, + 58.81167, + 11900, + 0, + 0, + 0, + 200e6, + -0.41, + 0, + 0, + 0, + 0, + ], # Cas ) ) @@ -80,11 +197,15 @@ def make_sky_model(sky0, settings, radius_deg, flux_min_outer_jy): sky_outer.num_sources, ) sky_outer.append(sky_inner) - LOG.info("Number of sources in output sky model: %d", sky_outer.num_sources) + LOG.info( + "Number of sources in output sky model: %d", sky_outer.num_sources + ) return sky_outer -def make_diff_image_stats(filename1, filename2, use_w_projection, out_image_root=None): +def make_diff_image_stats( + filename1, filename2, use_w_projection, out_image_root=None +): """Make an image of the difference between two visibility data sets. This function assumes that the observation parameters for both data sets @@ -98,7 +219,9 @@ def make_diff_image_stats(filename1, filename2, use_w_projection, out_image_root fov_ref_deg = 5.0 fov_deg = fov_ref_deg * (fov_ref_frequency_hz / frequency_hz) imager = oskar.Imager(precision="double") - imager.set(fov_deg=fov_deg, image_size=8192, fft_on_gpu=True, grid_on_gpu=True) + imager.set( + fov_deg=fov_deg, image_size=8192, fft_on_gpu=True, grid_on_gpu=True + ) if out_image_root is not None: imager.output_root = out_image_root @@ -144,10 +267,10 @@ def make_diff_image_stats(filename1, filename2, use_w_projection, out_image_root "image_medianabs": numpy.median(numpy.abs(image)), "image_mean": numpy.mean(image), "image_std": numpy.std(image), - "image_rms": numpy.sqrt(numpy.mean(image ** 2)), + "image_rms": numpy.sqrt(numpy.mean(image**2)), "image_centre_mean": numpy.mean(centre), "image_centre_std": numpy.std(centre), - "image_centre_rms": numpy.sqrt(numpy.mean(centre ** 2)), + "image_centre_rms": numpy.sqrt(numpy.mean(centre**2)), } @@ -155,7 +278,9 @@ def make_plot(prefix, field_name, metric_key, results, axis_freq): """Plot selected results.""" # Get data. data = numpy.zeros_like(axis_freq) - with numpy.nditer([axis_freq, data], op_flags=[["readonly"], ["writeonly"]]) as it: + with numpy.nditer( + [axis_freq, data], op_flags=[["readonly"], ["writeonly"]] + ) as it: for freq, d in it: key = "%s_%s_%d_MHz_iono" % (prefix, field_name, freq) if key in results: @@ -187,7 +312,9 @@ def run_single(prefix_field, settings, sky, freq_MHz, out0_name, results): LOG.info("Using cached results for '%s'", out) return settings["telescope/ionosphere_screen_type"] = "External" - settings["telescope/external_tec_screen/input_fits_file"] = "test_screen_60s.fits" + settings[ + "telescope/external_tec_screen/input_fits_file" + ] = "test_screen_60s.fits" settings["interferometer/oskar_vis_filename"] = out + ".vis" settings["interferometer/ms_filename"] = out + ".ms" make_vis_data(settings, sky) @@ -261,7 +388,12 @@ def run_set(prefix, base_settings, fields, axis_freq, plot_only): # Simulate the error case. run_single( - prefix_field, settings, sky, freq_MHz, out0 + ".vis", results + prefix_field, + settings, + sky, + freq_MHz, + out0 + ".vis", + results, ) # Generate plot for the field. diff --git a/scripts/run_sims_low_dd_effects.py b/scripts/run_sims_low_dd_effects.py index 41efa69ca6a32552ee13d610248f0c9e027f09fd..a1c9cebe155fc7ca61b4960920fe4fe2b40cd893 100644 --- a/scripts/run_sims_low_dd_effects.py +++ b/scripts/run_sims_low_dd_effects.py @@ -24,15 +24,132 @@ def bright_sources(): # Others from GLEAM reference paper, Hurley-Walker et al. (2017), Table 2. return numpy.array( ( - [50.67375, -37.20833, 528, 0, 0, 0, 178e6, -0.51, 0, 0, 0, 0], # For - [201.36667, -43.01917, 1370, 0, 0, 0, 200e6, -0.50, 0, 0, 0, 0], # Cen - [139.52500, -12.09556, 280, 0, 0, 0, 200e6, -0.96, 0, 0, 0, 0], # Hyd - [79.95833, -45.77889, 390, 0, 0, 0, 200e6, -0.99, 0, 0, 0, 0], # Pic - [252.78333, 4.99250, 377, 0, 0, 0, 200e6, -1.07, 0, 0, 0, 0], # Her - [187.70417, 12.39111, 861, 0, 0, 0, 200e6, -0.86, 0, 0, 0, 0], # Vir - [83.63333, 22.01444, 1340, 0, 0, 0, 200e6, -0.22, 0, 0, 0, 0], # Tau - [299.86667, 40.73389, 7920, 0, 0, 0, 200e6, -0.78, 0, 0, 0, 0], # Cyg - [350.86667, 58.81167, 11900, 0, 0, 0, 200e6, -0.41, 0, 0, 0, 0], # Cas + [ + 50.67375, + -37.20833, + 528, + 0, + 0, + 0, + 178e6, + -0.51, + 0, + 0, + 0, + 0, + ], # For + [ + 201.36667, + -43.01917, + 1370, + 0, + 0, + 0, + 200e6, + -0.50, + 0, + 0, + 0, + 0, + ], # Cen + [ + 139.52500, + -12.09556, + 280, + 0, + 0, + 0, + 200e6, + -0.96, + 0, + 0, + 0, + 0, + ], # Hyd + [ + 79.95833, + -45.77889, + 390, + 0, + 0, + 0, + 200e6, + -0.99, + 0, + 0, + 0, + 0, + ], # Pic + [ + 252.78333, + 4.99250, + 377, + 0, + 0, + 0, + 200e6, + -1.07, + 0, + 0, + 0, + 0, + ], # Her + [ + 187.70417, + 12.39111, + 861, + 0, + 0, + 0, + 200e6, + -0.86, + 0, + 0, + 0, + 0, + ], # Vir + [ + 83.63333, + 22.01444, + 1340, + 0, + 0, + 0, + 200e6, + -0.22, + 0, + 0, + 0, + 0, + ], # Tau + [ + 299.86667, + 40.73389, + 7920, + 0, + 0, + 0, + 200e6, + -0.78, + 0, + 0, + 0, + 0, + ], # Cyg + [ + 350.86667, + 58.81167, + 11900, + 0, + 0, + 0, + 200e6, + -0.41, + 0, + 0, + 0, + 0, + ], # Cas ) ) diff --git a/scripts/run_sims_obs_length.py b/scripts/run_sims_obs_length.py index e97b431b7bb1eea052bf1652a6b2bb668b459783..9bfea9c7b1605c11e484e0d33628f0aa3e50a083 100644 --- a/scripts/run_sims_obs_length.py +++ b/scripts/run_sims_obs_length.py @@ -104,11 +104,11 @@ def make_diff_image_stats(filename1, filename2, out_image_root=None): "image_maxabs": numpy.max(numpy.abs(image)), "image_medianabs": numpy.median(numpy.abs(image)), "image_mean": numpy.mean(image), - "image_rms": numpy.sqrt(numpy.mean(image ** 2)), + "image_rms": numpy.sqrt(numpy.mean(image**2)), "image_centre_maxabs": numpy.max(numpy.abs(centre)), "image_centre_mean": numpy.mean(centre), "image_centre_std": numpy.std(centre), - "image_centre_rms": numpy.sqrt(numpy.mean(centre ** 2)), + "image_centre_rms": numpy.sqrt(numpy.mean(centre**2)), } @@ -168,7 +168,9 @@ def make_sky_model(sky0, settings, radius_deg, flux_min_outer_jy): sky_outer.num_sources, ) sky_outer.append(sky_inner) - LOG.info("Number of sources in output sky model: %d", sky_outer.num_sources) + LOG.info( + "Number of sources in output sky model: %d", sky_outer.num_sources + ) return sky_outer @@ -279,7 +281,9 @@ def run_set( num_times = int(numpy.round(length / 60.0)) settings["observation/length"] = str(length) settings["observation/num_time_steps"] = str(num_times) - settings["observation/start_time_utc"] = get_start_time(ra_deg, length) + settings["observation/start_time_utc"] = get_start_time( + ra_deg, length + ) # Simulate the 'perfect' case. tel.override_element_gains(1.0, 0.0) @@ -357,7 +361,10 @@ def main(): "use_gpus": True, "max_sources_per_chunk": 23000, }, - "observation": {"start_frequency_hz": 100e6, "frequency_inc_hz": 100e3}, + "observation": { + "start_frequency_hz": 100e6, + "frequency_inc_hz": 100e3, + }, "telescope": { "input_directory": "SKA1-LOW_SKO-0000422_Rev3_38m.tm", "pol_mode": "Full", @@ -420,15 +427,132 @@ def main(): numpy.array( ( [266.41683, -29.00781, 2000, 0, 0, 0, 0, 0, 0, 3600, 3600, 0], - [50.67375, -37.20833, 528, 0, 0, 0, 178e6, -0.51, 0, 0, 0, 0], # For - [201.36667, -43.01917, 1370, 0, 0, 0, 200e6, -0.50, 0, 0, 0, 0], # Cen - [139.52500, -12.09556, 280, 0, 0, 0, 200e6, -0.96, 0, 0, 0, 0], # Hyd - [79.95833, -45.77889, 390, 0, 0, 0, 200e6, -0.99, 0, 0, 0, 0], # Pic - [252.78333, 4.99250, 377, 0, 0, 0, 200e6, -1.07, 0, 0, 0, 0], # Her - [187.70417, 12.39111, 861, 0, 0, 0, 200e6, -0.86, 0, 0, 0, 0], # Vir - [83.63333, 22.01444, 1340, 0, 0, 0, 200e6, -0.22, 0, 0, 0, 0], # Tau - [299.86667, 40.73389, 7920, 0, 0, 0, 200e6, -0.78, 0, 0, 0, 0], # Cyg - [350.86667, 58.81167, 11900, 0, 0, 0, 200e6, -0.41, 0, 0, 0, 0], # Cas + [ + 50.67375, + -37.20833, + 528, + 0, + 0, + 0, + 178e6, + -0.51, + 0, + 0, + 0, + 0, + ], # For + [ + 201.36667, + -43.01917, + 1370, + 0, + 0, + 0, + 200e6, + -0.50, + 0, + 0, + 0, + 0, + ], # Cen + [ + 139.52500, + -12.09556, + 280, + 0, + 0, + 0, + 200e6, + -0.96, + 0, + 0, + 0, + 0, + ], # Hyd + [ + 79.95833, + -45.77889, + 390, + 0, + 0, + 0, + 200e6, + -0.99, + 0, + 0, + 0, + 0, + ], # Pic + [ + 252.78333, + 4.99250, + 377, + 0, + 0, + 0, + 200e6, + -1.07, + 0, + 0, + 0, + 0, + ], # Her + [ + 187.70417, + 12.39111, + 861, + 0, + 0, + 0, + 200e6, + -0.86, + 0, + 0, + 0, + 0, + ], # Vir + [ + 83.63333, + 22.01444, + 1340, + 0, + 0, + 0, + 200e6, + -0.22, + 0, + 0, + 0, + 0, + ], # Tau + [ + 299.86667, + 40.73389, + 7920, + 0, + 0, + 0, + 200e6, + -0.78, + 0, + 0, + 0, + 0, + ], # Cyg + [ + 350.86667, + 58.81167, + 11900, + 0, + 0, + 0, + 200e6, + -0.41, + 0, + 0, + 0, + 0, + ], # Cas ) ) ) diff --git a/scripts/run_sims_spot_freqs_eda2.py b/scripts/run_sims_spot_freqs_eda2.py index c8765f70369d55a7d51d3e9cb78d75ce36ceae9f..f33dc0f53aa058ca0ddd1d3883f28ab6de77b5db 100644 --- a/scripts/run_sims_spot_freqs_eda2.py +++ b/scripts/run_sims_spot_freqs_eda2.py @@ -30,15 +30,132 @@ def bright_sources(): return numpy.array( ( [266.41683, -29.00781, 2000, 0, 0, 0, 0, 0, 0, 3600, 3600, 0], - [50.67375, -37.20833, 528, 0, 0, 0, 178e6, -0.51, 0, 0, 0, 0], # For - [201.36667, -43.01917, 1370, 0, 0, 0, 200e6, -0.50, 0, 0, 0, 0], # Cen - [139.52500, -12.09556, 280, 0, 0, 0, 200e6, -0.96, 0, 0, 0, 0], # Hyd - [79.95833, -45.77889, 390, 0, 0, 0, 200e6, -0.99, 0, 0, 0, 0], # Pic - [252.78333, 4.99250, 377, 0, 0, 0, 200e6, -1.07, 0, 0, 0, 0], # Her - [187.70417, 12.39111, 861, 0, 0, 0, 200e6, -0.86, 0, 0, 0, 0], # Vir - [83.63333, 22.01444, 1340, 0, 0, 0, 200e6, -0.22, 0, 0, 0, 0], # Tau - [299.86667, 40.73389, 7920, 0, 0, 0, 200e6, -0.78, 0, 0, 0, 0], # Cyg - [350.86667, 58.81167, 11900, 0, 0, 0, 200e6, -0.41, 0, 0, 0, 0], # Cas + [ + 50.67375, + -37.20833, + 528, + 0, + 0, + 0, + 178e6, + -0.51, + 0, + 0, + 0, + 0, + ], # For + [ + 201.36667, + -43.01917, + 1370, + 0, + 0, + 0, + 200e6, + -0.50, + 0, + 0, + 0, + 0, + ], # Cen + [ + 139.52500, + -12.09556, + 280, + 0, + 0, + 0, + 200e6, + -0.96, + 0, + 0, + 0, + 0, + ], # Hyd + [ + 79.95833, + -45.77889, + 390, + 0, + 0, + 0, + 200e6, + -0.99, + 0, + 0, + 0, + 0, + ], # Pic + [ + 252.78333, + 4.99250, + 377, + 0, + 0, + 0, + 200e6, + -1.07, + 0, + 0, + 0, + 0, + ], # Her + [ + 187.70417, + 12.39111, + 861, + 0, + 0, + 0, + 200e6, + -0.86, + 0, + 0, + 0, + 0, + ], # Vir + [ + 83.63333, + 22.01444, + 1340, + 0, + 0, + 0, + 200e6, + -0.22, + 0, + 0, + 0, + 0, + ], # Tau + [ + 299.86667, + 40.73389, + 7920, + 0, + 0, + 0, + 200e6, + -0.78, + 0, + 0, + 0, + 0, + ], # Cyg + [ + 350.86667, + 58.81167, + 11900, + 0, + 0, + 0, + 200e6, + -0.41, + 0, + 0, + 0, + 0, + ], # Cas ) ) @@ -80,11 +197,15 @@ def make_sky_model(sky0, settings, radius_deg, flux_min_outer_jy): sky_outer.num_sources, ) sky_outer.append(sky_inner) - LOG.info("Number of sources in output sky model: %d", sky_outer.num_sources) + LOG.info( + "Number of sources in output sky model: %d", sky_outer.num_sources + ) return sky_outer -def make_diff_image_stats(filename1, filename2, use_w_projection, out_image_root=None): +def make_diff_image_stats( + filename1, filename2, use_w_projection, out_image_root=None +): """Make an image of the difference between two visibility data sets. This function assumes that the observation parameters for both data sets @@ -98,7 +219,9 @@ def make_diff_image_stats(filename1, filename2, use_w_projection, out_image_root fov_ref_deg = 5.0 fov_deg = fov_ref_deg * (fov_ref_frequency_hz / frequency_hz) imager = oskar.Imager(precision="double") - imager.set(fov_deg=fov_deg, image_size=8192, fft_on_gpu=True, grid_on_gpu=True) + imager.set( + fov_deg=fov_deg, image_size=8192, fft_on_gpu=True, grid_on_gpu=True + ) if out_image_root is not None: imager.output_root = out_image_root @@ -144,10 +267,10 @@ def make_diff_image_stats(filename1, filename2, use_w_projection, out_image_root "image_medianabs": numpy.median(numpy.abs(image)), "image_mean": numpy.mean(image), "image_std": numpy.std(image), - "image_rms": numpy.sqrt(numpy.mean(image ** 2)), + "image_rms": numpy.sqrt(numpy.mean(image**2)), "image_centre_mean": numpy.mean(centre), "image_centre_std": numpy.std(centre), - "image_centre_rms": numpy.sqrt(numpy.mean(centre ** 2)), + "image_centre_rms": numpy.sqrt(numpy.mean(centre**2)), } @@ -265,7 +388,9 @@ def run_set(prefix, base_settings, fields, axis_freq, axis_gain, plot_only): ra_deg = float(settings["observation/phase_centre_ra_deg"]) dec_deg = float(settings["observation/phase_centre_dec_deg"]) length_sec = float(settings["observation/length"]) - settings["observation/start_frequency_hz"] = str(freq_MHz * 1e6) + settings["observation/start_frequency_hz"] = str( + freq_MHz * 1e6 + ) settings["observation/start_time_utc"] = get_start_time( ra_deg, length_sec ) @@ -280,7 +405,10 @@ def run_set(prefix, base_settings, fields, axis_freq, axis_gain, plot_only): # Simulate the 'perfect' case. tel.override_element_gains(1.0, 0.0) tel.override_element_cable_length_errors(0.0) - out0_name = "%s_%d_MHz_no_errors.vis" % (prefix_field, freq_MHz) + out0_name = "%s_%d_MHz_no_errors.vis" % ( + prefix_field, + freq_MHz, + ) settings["interferometer/oskar_vis_filename"] = out0_name make_vis_data(settings, sky, tel) @@ -298,8 +426,22 @@ def run_set(prefix, base_settings, fields, axis_freq, axis_gain, plot_only): ) # Generate plot for the field. - make_plot(prefix, field_name, "image_centre_rms", results, axis_freq, axis_gain) - make_plot(prefix, field_name, "image_medianabs", results, axis_freq, axis_gain) + make_plot( + prefix, + field_name, + "image_centre_rms", + results, + axis_freq, + axis_gain, + ) + make_plot( + prefix, + field_name, + "image_medianabs", + results, + axis_freq, + axis_gain, + ) # Save result set. with open(json_file, "w") as output_file: @@ -358,7 +500,9 @@ def main(): # GLEAM + A-team sky model simulations. plot_only = True - run_set("GLEAM_A-team", base_settings, fields, axis_freq, axis_gain, plot_only) + run_set( + "GLEAM_A-team", base_settings, fields, axis_freq, axis_gain, plot_only + ) if __name__ == "__main__": diff --git a/scripts/run_sims_spot_freqs_skala4.py b/scripts/run_sims_spot_freqs_skala4.py index 717d83cd6d7a5ebc92ad0320b349afbcc8091817..96ed829b72535a8a7391c31e822cfc112ce5bd64 100644 --- a/scripts/run_sims_spot_freqs_skala4.py +++ b/scripts/run_sims_spot_freqs_skala4.py @@ -30,15 +30,132 @@ def bright_sources(): return numpy.array( ( [266.41683, -29.00781, 2000, 0, 0, 0, 0, 0, 0, 3600, 3600, 0], - [50.67375, -37.20833, 528, 0, 0, 0, 178e6, -0.51, 0, 0, 0, 0], # For - [201.36667, -43.01917, 1370, 0, 0, 0, 200e6, -0.50, 0, 0, 0, 0], # Cen - [139.52500, -12.09556, 280, 0, 0, 0, 200e6, -0.96, 0, 0, 0, 0], # Hyd - [79.95833, -45.77889, 390, 0, 0, 0, 200e6, -0.99, 0, 0, 0, 0], # Pic - [252.78333, 4.99250, 377, 0, 0, 0, 200e6, -1.07, 0, 0, 0, 0], # Her - [187.70417, 12.39111, 861, 0, 0, 0, 200e6, -0.86, 0, 0, 0, 0], # Vir - [83.63333, 22.01444, 1340, 0, 0, 0, 200e6, -0.22, 0, 0, 0, 0], # Tau - [299.86667, 40.73389, 7920, 0, 0, 0, 200e6, -0.78, 0, 0, 0, 0], # Cyg - [350.86667, 58.81167, 11900, 0, 0, 0, 200e6, -0.41, 0, 0, 0, 0], # Cas + [ + 50.67375, + -37.20833, + 528, + 0, + 0, + 0, + 178e6, + -0.51, + 0, + 0, + 0, + 0, + ], # For + [ + 201.36667, + -43.01917, + 1370, + 0, + 0, + 0, + 200e6, + -0.50, + 0, + 0, + 0, + 0, + ], # Cen + [ + 139.52500, + -12.09556, + 280, + 0, + 0, + 0, + 200e6, + -0.96, + 0, + 0, + 0, + 0, + ], # Hyd + [ + 79.95833, + -45.77889, + 390, + 0, + 0, + 0, + 200e6, + -0.99, + 0, + 0, + 0, + 0, + ], # Pic + [ + 252.78333, + 4.99250, + 377, + 0, + 0, + 0, + 200e6, + -1.07, + 0, + 0, + 0, + 0, + ], # Her + [ + 187.70417, + 12.39111, + 861, + 0, + 0, + 0, + 200e6, + -0.86, + 0, + 0, + 0, + 0, + ], # Vir + [ + 83.63333, + 22.01444, + 1340, + 0, + 0, + 0, + 200e6, + -0.22, + 0, + 0, + 0, + 0, + ], # Tau + [ + 299.86667, + 40.73389, + 7920, + 0, + 0, + 0, + 200e6, + -0.78, + 0, + 0, + 0, + 0, + ], # Cyg + [ + 350.86667, + 58.81167, + 11900, + 0, + 0, + 0, + 200e6, + -0.41, + 0, + 0, + 0, + 0, + ], # Cas ) ) @@ -80,11 +197,15 @@ def make_sky_model(sky0, settings, radius_deg, flux_min_outer_jy): sky_outer.num_sources, ) sky_outer.append(sky_inner) - LOG.info("Number of sources in output sky model: %d", sky_outer.num_sources) + LOG.info( + "Number of sources in output sky model: %d", sky_outer.num_sources + ) return sky_outer -def make_diff_image_stats(filename1, filename2, use_w_projection, out_image_root=None): +def make_diff_image_stats( + filename1, filename2, use_w_projection, out_image_root=None +): """Make an image of the difference between two visibility data sets. This function assumes that the observation parameters for both data sets @@ -98,7 +219,9 @@ def make_diff_image_stats(filename1, filename2, use_w_projection, out_image_root fov_ref_deg = 5.0 fov_deg = fov_ref_deg * (fov_ref_frequency_hz / frequency_hz) imager = oskar.Imager(precision="double") - imager.set(fov_deg=fov_deg, image_size=8192, fft_on_gpu=True, grid_on_gpu=True) + imager.set( + fov_deg=fov_deg, image_size=8192, fft_on_gpu=True, grid_on_gpu=True + ) if out_image_root is not None: imager.output_root = out_image_root @@ -144,10 +267,10 @@ def make_diff_image_stats(filename1, filename2, use_w_projection, out_image_root "image_medianabs": numpy.median(numpy.abs(image)), "image_mean": numpy.mean(image), "image_std": numpy.std(image), - "image_rms": numpy.sqrt(numpy.mean(image ** 2)), + "image_rms": numpy.sqrt(numpy.mean(image**2)), "image_centre_mean": numpy.mean(centre), "image_centre_std": numpy.std(centre), - "image_centre_rms": numpy.sqrt(numpy.mean(centre ** 2)), + "image_centre_rms": numpy.sqrt(numpy.mean(centre**2)), } @@ -265,7 +388,9 @@ def run_set(prefix, base_settings, fields, axis_freq, axis_gain, plot_only): ra_deg = float(settings["observation/phase_centre_ra_deg"]) dec_deg = float(settings["observation/phase_centre_dec_deg"]) length_sec = float(settings["observation/length"]) - settings["observation/start_frequency_hz"] = str(freq_MHz * 1e6) + settings["observation/start_frequency_hz"] = str( + freq_MHz * 1e6 + ) settings["observation/start_time_utc"] = get_start_time( ra_deg, length_sec ) @@ -280,7 +405,10 @@ def run_set(prefix, base_settings, fields, axis_freq, axis_gain, plot_only): # Simulate the 'perfect' case. tel.override_element_gains(1.0, 0.0) tel.override_element_cable_length_errors(0.0) - out0_name = "%s_%d_MHz_no_errors.vis" % (prefix_field, freq_MHz) + out0_name = "%s_%d_MHz_no_errors.vis" % ( + prefix_field, + freq_MHz, + ) settings["interferometer/oskar_vis_filename"] = out0_name make_vis_data(settings, sky, tel) @@ -298,8 +426,22 @@ def run_set(prefix, base_settings, fields, axis_freq, axis_gain, plot_only): ) # Generate plot for the field. - make_plot(prefix, field_name, "image_centre_rms", results, axis_freq, axis_gain) - make_plot(prefix, field_name, "image_medianabs", results, axis_freq, axis_gain) + make_plot( + prefix, + field_name, + "image_centre_rms", + results, + axis_freq, + axis_gain, + ) + make_plot( + prefix, + field_name, + "image_medianabs", + results, + axis_freq, + axis_gain, + ) # Save result set. with open(json_file, "w") as output_file: @@ -358,7 +500,9 @@ def main(): # GLEAM + A-team sky model simulations. plot_only = True - run_set("GLEAM_A-team", base_settings, fields, axis_freq, axis_gain, plot_only) + run_set( + "GLEAM_A-team", base_settings, fields, axis_freq, axis_gain, plot_only + ) if __name__ == "__main__": diff --git a/scripts/run_sims_station_layouts.py b/scripts/run_sims_station_layouts.py index 9d97204ede3f8ec79a3b74d71cc2d1a582809fd4..49267fd1f6e823d535eec9393fbc55a7e0566ce4 100644 --- a/scripts/run_sims_station_layouts.py +++ b/scripts/run_sims_station_layouts.py @@ -31,15 +31,132 @@ def bright_sources(): return numpy.array( ( [266.41683, -29.00781, 2000, 0, 0, 0, 0, 0, 0, 3600, 3600, 0], - [50.67375, -37.20833, 528, 0, 0, 0, 178e6, -0.51, 0, 0, 0, 0], # For - [201.36667, -43.01917, 1370, 0, 0, 0, 200e6, -0.50, 0, 0, 0, 0], # Cen - [139.52500, -12.09556, 280, 0, 0, 0, 200e6, -0.96, 0, 0, 0, 0], # Hyd - [79.95833, -45.77889, 390, 0, 0, 0, 200e6, -0.99, 0, 0, 0, 0], # Pic - [252.78333, 4.99250, 377, 0, 0, 0, 200e6, -1.07, 0, 0, 0, 0], # Her - [187.70417, 12.39111, 861, 0, 0, 0, 200e6, -0.86, 0, 0, 0, 0], # Vir - [83.63333, 22.01444, 1340, 0, 0, 0, 200e6, -0.22, 0, 0, 0, 0], # Tau - [299.86667, 40.73389, 7920, 0, 0, 0, 200e6, -0.78, 0, 0, 0, 0], # Cyg - [350.86667, 58.81167, 11900, 0, 0, 0, 200e6, -0.41, 0, 0, 0, 0], # Cas + [ + 50.67375, + -37.20833, + 528, + 0, + 0, + 0, + 178e6, + -0.51, + 0, + 0, + 0, + 0, + ], # For + [ + 201.36667, + -43.01917, + 1370, + 0, + 0, + 0, + 200e6, + -0.50, + 0, + 0, + 0, + 0, + ], # Cen + [ + 139.52500, + -12.09556, + 280, + 0, + 0, + 0, + 200e6, + -0.96, + 0, + 0, + 0, + 0, + ], # Hyd + [ + 79.95833, + -45.77889, + 390, + 0, + 0, + 0, + 200e6, + -0.99, + 0, + 0, + 0, + 0, + ], # Pic + [ + 252.78333, + 4.99250, + 377, + 0, + 0, + 0, + 200e6, + -1.07, + 0, + 0, + 0, + 0, + ], # Her + [ + 187.70417, + 12.39111, + 861, + 0, + 0, + 0, + 200e6, + -0.86, + 0, + 0, + 0, + 0, + ], # Vir + [ + 83.63333, + 22.01444, + 1340, + 0, + 0, + 0, + 200e6, + -0.22, + 0, + 0, + 0, + 0, + ], # Tau + [ + 299.86667, + 40.73389, + 7920, + 0, + 0, + 0, + 200e6, + -0.78, + 0, + 0, + 0, + 0, + ], # Cyg + [ + 350.86667, + 58.81167, + 11900, + 0, + 0, + 0, + 200e6, + -0.41, + 0, + 0, + 0, + 0, + ], # Cas ) ) @@ -80,11 +197,15 @@ def make_sky_model(sky0, settings, radius_deg, flux_min_outer_jy): sky_outer.num_sources, ) sky_outer.append(sky_inner) - LOG.info("Number of sources in output sky model: %d", sky_outer.num_sources) + LOG.info( + "Number of sources in output sky model: %d", sky_outer.num_sources + ) return sky_outer -def make_diff_image_stats(filename1, filename2, use_w_projection, out_image_root=None): +def make_diff_image_stats( + filename1, filename2, use_w_projection, out_image_root=None +): """Make an image of the difference between two visibility data sets. This function assumes that the observation parameters for both data sets @@ -98,7 +219,9 @@ def make_diff_image_stats(filename1, filename2, use_w_projection, out_image_root fov_ref_deg = 5.0 fov_deg = fov_ref_deg * (fov_ref_frequency_hz / frequency_hz) imager = oskar.Imager(precision="double") - imager.set(fov_deg=fov_deg, image_size=8192, fft_on_gpu=True, grid_on_gpu=True) + imager.set( + fov_deg=fov_deg, image_size=8192, fft_on_gpu=True, grid_on_gpu=True + ) if out_image_root is not None: imager.output_root = out_image_root @@ -144,10 +267,10 @@ def make_diff_image_stats(filename1, filename2, use_w_projection, out_image_root "image_medianabs": numpy.median(numpy.abs(image)), "image_mean": numpy.mean(image), "image_std": numpy.std(image), - "image_rms": numpy.sqrt(numpy.mean(image ** 2)), + "image_rms": numpy.sqrt(numpy.mean(image**2)), "image_centre_mean": numpy.mean(centre), "image_centre_std": numpy.std(centre), - "image_centre_rms": numpy.sqrt(numpy.mean(centre ** 2)), + "image_centre_rms": numpy.sqrt(numpy.mean(centre**2)), } @@ -188,7 +311,9 @@ def make_plots(prefix, metric_key, results, fields, axis_stations, axis_freq): plt.close("all") -def run_single(prefix_field, settings, sky, stations, freq_MHz, out0_name, results): +def run_single( + prefix_field, settings, sky, stations, freq_MHz, out0_name, results +): """Run a single simulation and generate image statistics for it.""" out = "%s_%03d_stations_%d_MHz" % (prefix_field, stations, freq_MHz) if out in results: @@ -241,7 +366,9 @@ def run_set(prefix, base_settings, fields, axis_stations, axis_freq, results): ra_deg = float(settings["observation/phase_centre_ra_deg"]) length_sec = float(settings["observation/length"]) settings["observation/start_frequency_hz"] = str(freq_MHz * 1e6) - settings["observation/start_time_utc"] = get_start_time(ra_deg, length_sec) + settings["observation/start_time_utc"] = get_start_time( + ra_deg, length_sec + ) # Create the sky model. sky = make_sky_model(sky0, settings, 20.0, 10.0) @@ -253,7 +380,10 @@ def run_set(prefix, base_settings, fields, axis_stations, axis_freq, results): settings[ "telescope/input_directory" ] = "SKA1-LOW_SKO-0000422_Rev3_SKALA4_REF_512_different_38m_stations.tm" - out0_name = "%s_REF_512_stations_%d_MHz.vis" % (prefix_field, freq_MHz) + out0_name = "%s_REF_512_stations_%d_MHz.vis" % ( + prefix_field, + freq_MHz, + ) settings["interferometer/oskar_vis_filename"] = out0_name make_vis_data(settings, sky) @@ -264,7 +394,13 @@ def run_set(prefix, base_settings, fields, axis_stations, axis_freq, results): % stations ) run_single( - prefix_field, settings, sky, stations, freq_MHz, out0_name, results + prefix_field, + settings, + sky, + stations, + freq_MHz, + out0_name, + results, ) # Update saved result set. @@ -319,8 +455,12 @@ def main(): run_set(prefix, base_settings, fields, axis_stations, axis_freq, results) # Generate plots. - make_plots(prefix, "image_centre_rms", results, fields, axis_stations, axis_freq) - make_plots(prefix, "image_medianabs", results, fields, axis_stations, axis_freq) + make_plots( + prefix, "image_centre_rms", results, fields, axis_stations, axis_freq + ) + make_plots( + prefix, "image_medianabs", results, fields, axis_stations, axis_freq + ) if __name__ == "__main__": diff --git a/scripts/run_sims_station_layouts_bp.py b/scripts/run_sims_station_layouts_bp.py index 782b848a6c146fb01ca5abcb7c77f842426e8b7a..5862a2a47e56622b91db94c2f032b844d64e752f 100644 --- a/scripts/run_sims_station_layouts_bp.py +++ b/scripts/run_sims_station_layouts_bp.py @@ -32,7 +32,9 @@ def run_set(prefix, base_settings, fields, axis_stations, axis_freq): settings.from_dict(settings_dict) ra_deg = float(settings["observation/phase_centre_ra_deg"]) length_sec = float(settings["observation/length"]) - settings["observation/start_time_utc"] = get_start_time(ra_deg, length_sec) + settings["observation/start_time_utc"] = get_start_time( + ra_deg, length_sec + ) # Iterate over frequencies. for freq_MHz in axis_freq: @@ -51,7 +53,10 @@ def run_set(prefix, base_settings, fields, axis_stations, axis_freq): freq_MHz, ) test_path = "" - if settings["beam_pattern/output/separate_time_and_channel"] == "True": + if ( + settings["beam_pattern/output/separate_time_and_channel"] + == "True" + ): test_path = root_path + "*TIME_SEP*" else: test_path = root_path + "*TIME_AVG*" @@ -109,7 +114,13 @@ def main(): axis_freq = [110, 230] # axis_stations = [1, 4, 16, 512] axis_stations = [1] - run_set("auto_power_beam_snapshot", base_settings, fields, axis_stations, axis_freq) + run_set( + "auto_power_beam_snapshot", + base_settings, + fields, + axis_stations, + axis_freq, + ) if __name__ == "__main__": diff --git a/scripts/run_sims_station_layouts_fscn.py b/scripts/run_sims_station_layouts_fscn.py index 6352375e5df600daa2ca001ee0e654b29a6fe3dd..29d7387dc4b6b6de2d1fed10c6d9b55836306ca9 100644 --- a/scripts/run_sims_station_layouts_fscn.py +++ b/scripts/run_sims_station_layouts_fscn.py @@ -36,15 +36,132 @@ def bright_sources(): return numpy.array( ( [266.41683, -29.00781, 2000, 0, 0, 0, 0, 0, 0, 3600, 3600, 0], - [50.67375, -37.20833, 528, 0, 0, 0, 178e6, -0.51, 0, 0, 0, 0], # For - [201.36667, -43.01917, 1370, 0, 0, 0, 200e6, -0.50, 0, 0, 0, 0], # Cen - [139.52500, -12.09556, 280, 0, 0, 0, 200e6, -0.96, 0, 0, 0, 0], # Hyd - [79.95833, -45.77889, 390, 0, 0, 0, 200e6, -0.99, 0, 0, 0, 0], # Pic - [252.78333, 4.99250, 377, 0, 0, 0, 200e6, -1.07, 0, 0, 0, 0], # Her - [187.70417, 12.39111, 861, 0, 0, 0, 200e6, -0.86, 0, 0, 0, 0], # Vir - [83.63333, 22.01444, 1340, 0, 0, 0, 200e6, -0.22, 0, 0, 0, 0], # Tau - [299.86667, 40.73389, 7920, 0, 0, 0, 200e6, -0.78, 0, 0, 0, 0], # Cyg - [350.86667, 58.81167, 11900, 0, 0, 0, 200e6, -0.41, 0, 0, 0, 0], # Cas + [ + 50.67375, + -37.20833, + 528, + 0, + 0, + 0, + 178e6, + -0.51, + 0, + 0, + 0, + 0, + ], # For + [ + 201.36667, + -43.01917, + 1370, + 0, + 0, + 0, + 200e6, + -0.50, + 0, + 0, + 0, + 0, + ], # Cen + [ + 139.52500, + -12.09556, + 280, + 0, + 0, + 0, + 200e6, + -0.96, + 0, + 0, + 0, + 0, + ], # Hyd + [ + 79.95833, + -45.77889, + 390, + 0, + 0, + 0, + 200e6, + -0.99, + 0, + 0, + 0, + 0, + ], # Pic + [ + 252.78333, + 4.99250, + 377, + 0, + 0, + 0, + 200e6, + -1.07, + 0, + 0, + 0, + 0, + ], # Her + [ + 187.70417, + 12.39111, + 861, + 0, + 0, + 0, + 200e6, + -0.86, + 0, + 0, + 0, + 0, + ], # Vir + [ + 83.63333, + 22.01444, + 1340, + 0, + 0, + 0, + 200e6, + -0.22, + 0, + 0, + 0, + 0, + ], # Tau + [ + 299.86667, + 40.73389, + 7920, + 0, + 0, + 0, + 200e6, + -0.78, + 0, + 0, + 0, + 0, + ], # Cyg + [ + 350.86667, + 58.81167, + 11900, + 0, + 0, + 0, + 200e6, + -0.41, + 0, + 0, + 0, + 0, + ], # Cas ) ) @@ -121,10 +238,10 @@ def make_image_stats(filename, use_w_projection, out_image_root=None): "image_medianabs": numpy.median(numpy.abs(image)), "image_mean": numpy.mean(image), "image_std": numpy.std(image), - "image_rms": numpy.sqrt(numpy.mean(image ** 2)), + "image_rms": numpy.sqrt(numpy.mean(image**2)), "image_centre_mean": numpy.mean(centre), "image_centre_std": numpy.std(centre), - "image_centre_rms": numpy.sqrt(numpy.mean(centre ** 2)), + "image_centre_rms": numpy.sqrt(numpy.mean(centre**2)), } @@ -214,7 +331,9 @@ def run_set(prefix, base_settings, fields, axis_stations, axis_freq, results): settings.from_dict(settings_dict) ra_deg = float(settings["observation/phase_centre_ra_deg"]) length_sec = float(settings["observation/length"]) - settings["observation/start_time_utc"] = get_start_time(ra_deg, length_sec) + settings["observation/start_time_utc"] = get_start_time( + ra_deg, length_sec + ) # Iterate over frequencies. for freq_MHz in axis_freq: @@ -235,7 +354,9 @@ def run_set(prefix, base_settings, fields, axis_stations, axis_freq, results): % stations ) t0 = time.time() - run_single(prefix_field, settings, sky, stations, freq_MHz, results) + run_single( + prefix_field, settings, sky, stations, freq_MHz, results + ) duration = time.time() - t0 # Update saved result set. @@ -309,8 +430,12 @@ def main(): run_set(prefix, base_settings, fields, axis_stations, axis_freq, results) # Generate final plots. - make_plots(prefix, "image_centre_rms", results, fields, axis_stations, axis_freq) - make_plots(prefix, "image_medianabs", results, fields, axis_stations, axis_freq) + make_plots( + prefix, "image_centre_rms", results, fields, axis_stations, axis_freq + ) + make_plots( + prefix, "image_medianabs", results, fields, axis_stations, axis_freq + ) if __name__ == "__main__": diff --git a/scripts/run_sims_unpolarised_grid_beam_errors.py b/scripts/run_sims_unpolarised_grid_beam_errors.py index 883e805291824ee17c00c12aa0e00709e5afac18..a69e8a6fba3f5eba25400d81e77a94d3e4489edc 100644 --- a/scripts/run_sims_unpolarised_grid_beam_errors.py +++ b/scripts/run_sims_unpolarised_grid_beam_errors.py @@ -28,7 +28,9 @@ class Simulator(oskar.Interferometer): method. """ - def __init__(self, imager=None, precision=None, settings=None, ref_vis=None): + def __init__( + self, imager=None, precision=None, settings=None, ref_vis=None + ): """Creates the simulator, storing a handle to the imager. Args: @@ -69,10 +71,14 @@ class Simulator(oskar.Interferometer): if not self.coords_only: if self._ref_vis: # Read the reference visibility data block. - self._ref_block.read(self._ref_hdr, self._ref_handle, block_index) + self._ref_block.read( + self._ref_hdr, self._ref_handle, block_index + ) # Subtract reference block from this one. - block.cross_correlations()[:] -= self._ref_block.cross_correlations()[:] + block.cross_correlations()[ + : + ] -= self._ref_block.cross_correlations()[:] # Write (modified) block to any open files. self.write_block(block, block_index) @@ -138,7 +144,12 @@ def make_plots(prefix, ra0, dec0, freq_MHz, sky_pos, results): # Loop over source positions. for s in range(num_pos): - key = "%s_%.1f_%.1f_%.0f_MHz" % (prefix, sky_pos[s, 0], sky_pos[s, 1], freq_MHz) + key = "%s_%.1f_%.1f_%.0f_MHz" % ( + prefix, + sky_pos[s, 0], + sky_pos[s, 1], + freq_MHz, + ) if key in results: ra = numpy.append(ra, numpy.radians(sky_pos[s, 0])) dec = numpy.append(dec, numpy.radians(sky_pos[s, 1])) @@ -166,7 +177,9 @@ def make_plots(prefix, ra0, dec0, freq_MHz, sky_pos, results): # Common axis labels. ax = fig.add_subplot(111, frameon=False) - plt.tick_params(labelcolor="none", top="off", bottom="off", left="off", right="off") + plt.tick_params( + labelcolor="none", top="off", bottom="off", left="off", right="off" + ) ax.set_xlabel("x direction cosine", labelpad=20) ax.set_ylabel("y direction cosine", labelpad=40) # Hide ticks for this new axis. @@ -179,7 +192,9 @@ def make_plots(prefix, ra0, dec0, freq_MHz, sky_pos, results): plt.close("all") -def run_set(prefix, axis_freq, ra0_deg, dec0_deg, base_settings, sky_pos, results): +def run_set( + prefix, axis_freq, ra0_deg, dec0_deg, base_settings, sky_pos, results +): """Run a set of simulations.""" settings = oskar.SettingsTree("oskar_sim_interferometer") settings.from_dict(base_settings) @@ -229,7 +244,11 @@ def run_set(prefix, axis_freq, ra0_deg, dec0_deg, base_settings, sky_pos, result # Run simulation and make images of all four polarisations. settings["interferometer/oskar_vis_filename"] = "" - LOG.info("Running simulation for source at (%.3f, %.3f)", ra_deg, dec_deg) + LOG.info( + "Running simulation for source at (%.3f, %.3f)", + ra_deg, + dec_deg, + ) sim = Simulator(imager, settings=settings, ref_vis=ref_vis) sim.set_sky_model(sky) sim.set_telescope_model(tel_errors) @@ -298,7 +317,9 @@ def main(): axis_freq = [50, 70, 110, 137, 160, 230, 320] # axis_freq = [230] prefix = "pol_scan_beam_errors_stokes_i_skala4" - run_set(prefix, axis_freq, ra0_deg, dec0_deg, base_settings, sky_pos, results) + run_set( + prefix, axis_freq, ra0_deg, dec0_deg, base_settings, sky_pos, results + ) if __name__ == "__main__": diff --git a/scripts/run_sims_unpolarised_grid_element.py b/scripts/run_sims_unpolarised_grid_element.py index 39e8709c9e4553003627c2ab79c5343604c6a69f..5bad2fa7c5abf1f76ca7a3915e7d6d45259b743f 100644 --- a/scripts/run_sims_unpolarised_grid_element.py +++ b/scripts/run_sims_unpolarised_grid_element.py @@ -150,7 +150,9 @@ def make_plots(prefix, ra0, dec0, frequency_hz, sky_pos, results): # Common axis labels. ax = fig.add_subplot(111, frameon=False) - plt.tick_params(labelcolor="none", top="off", bottom="off", left="off", right="off") + plt.tick_params( + labelcolor="none", top="off", bottom="off", left="off", right="off" + ) ax.set_xlabel("x direction cosine", labelpad=20) ax.set_ylabel("y direction cosine", labelpad=40) # Hide ticks for this new axis. @@ -174,7 +176,12 @@ def run_set(prefix, ra0_deg, dec0_deg, base_settings, sky_pos, results): dec_deg = sky_pos[s, 1] # Check if already done. - key = "%s_%.1f_%.1f_%.0f_MHz" % (prefix, ra_deg, dec_deg, frequency_hz * 1e-6) + key = "%s_%.1f_%.1f_%.0f_MHz" % ( + prefix, + ra_deg, + dec_deg, + frequency_hz * 1e-6, + ) if key in results: continue @@ -192,7 +199,9 @@ def run_set(prefix, ra0_deg, dec0_deg, base_settings, sky_pos, results): imager.output_root = key # Run simulation and make images of all four polarisations. - LOG.info("Running simulation for source at (%.3f, %.3f)", ra_deg, dec_deg) + LOG.info( + "Running simulation for source at (%.3f, %.3f)", ra_deg, dec_deg + ) sim = Simulator(imager, settings=settings) sim.set_sky_model(sky) output = sim.run(return_images=4) diff --git a/scripts/run_sims_unpolarised_grid_full.py b/scripts/run_sims_unpolarised_grid_full.py index c016e8acf166ab4ad3b6f6c6da7a200fe0401989..306a0320aac45b939885cae7f1de5e443f31ee21 100644 --- a/scripts/run_sims_unpolarised_grid_full.py +++ b/scripts/run_sims_unpolarised_grid_full.py @@ -163,7 +163,9 @@ def make_plots(prefix, ra0, dec0, frequency_hz, sky_pos, results): # Common axis labels. ax = fig.add_subplot(111, frameon=False) - plt.tick_params(labelcolor="none", top="off", bottom="off", left="off", right="off") + plt.tick_params( + labelcolor="none", top="off", bottom="off", left="off", right="off" + ) ax.set_xlabel("x direction cosine", labelpad=20) ax.set_ylabel("y direction cosine", labelpad=40) # Hide ticks for this new axis. @@ -179,7 +181,9 @@ def make_plots(prefix, ra0, dec0, frequency_hz, sky_pos, results): # normalised I, pol angle, fractional linear, fractional circular. norm_i = stokes_i / numpy.max(stokes_i) pol_angle = (0.5 * numpy.arctan2(stokes_u, stokes_q)) * 180.0 / numpy.pi - frac_linear = numpy.sqrt(numpy.square(stokes_q) + numpy.square(stokes_u)) / stokes_i + frac_linear = ( + numpy.sqrt(numpy.square(stokes_q) + numpy.square(stokes_u)) / stokes_i + ) frac_circular = numpy.abs(stokes_v) / stokes_i fig = plt.figure(figsize=(8, 6)) @@ -188,13 +192,17 @@ def make_plots(prefix, ra0, dec0, frequency_hz, sky_pos, results): ax = fig.add_subplot(222) plot_panel(ax, x, y, frac_linear, "Measured Linear Polarisation Fraction") ax = fig.add_subplot(223) - plot_panel(ax, x, y, frac_circular, "Measured Circular Polarisation Fraction") + plot_panel( + ax, x, y, frac_circular, "Measured Circular Polarisation Fraction" + ) ax = fig.add_subplot(224) plot_panel(ax, x, y, pol_angle, "Measured Polarisation Angle [deg]") # Common axis labels. ax = fig.add_subplot(111, frameon=False) - plt.tick_params(labelcolor="none", top="off", bottom="off", left="off", right="off") + plt.tick_params( + labelcolor="none", top="off", bottom="off", left="off", right="off" + ) ax.set_xlabel("x direction cosine", labelpad=20) ax.set_ylabel("y direction cosine", labelpad=40) # Hide ticks for this new axis. @@ -207,7 +215,9 @@ def make_plots(prefix, ra0, dec0, frequency_hz, sky_pos, results): plt.close("all") -def run_set(prefix, axis_freq, ra0_deg, dec0_deg, base_settings, sky_pos, results): +def run_set( + prefix, axis_freq, ra0_deg, dec0_deg, base_settings, sky_pos, results +): """Run a set of simulations.""" settings = oskar.SettingsTree("oskar_sim_interferometer") num_pos = sky_pos.shape[0] @@ -241,7 +251,11 @@ def run_set(prefix, axis_freq, ra0_deg, dec0_deg, base_settings, sky_pos, result # Run simulation and make images of all four polarisations. settings["observation/phase_centre_ra_deg"] = str(ra_deg) settings["observation/phase_centre_dec_deg"] = str(dec_deg) - LOG.info("Running simulation for source at (%.3f, %.3f)", ra_deg, dec_deg) + LOG.info( + "Running simulation for source at (%.3f, %.3f)", + ra_deg, + dec_deg, + ) sim = Simulator(imager, settings=settings) sim.set_sky_model(sky) output = sim.run(return_images=4) @@ -312,7 +326,9 @@ def main(): # axis_freq = [50, 70, 110, 137, 160, 230, 320] axis_freq = [230] prefix = "pol_scan_stokes_i_dipole_no_parangle" - run_set(prefix, axis_freq, ra0_deg, dec0_deg, base_settings, sky_pos, results) + run_set( + prefix, axis_freq, ra0_deg, dec0_deg, base_settings, sky_pos, results + ) if __name__ == "__main__": diff --git a/scripts/run_sims_wide_bandwidth.py b/scripts/run_sims_wide_bandwidth.py index 4d53e9cf1b17c2ac1513fddd7141f298b237f7a8..26b1e49b1e4894c367366e3d153a5f9259a219ec 100644 --- a/scripts/run_sims_wide_bandwidth.py +++ b/scripts/run_sims_wide_bandwidth.py @@ -31,15 +31,132 @@ def bright_sources(): return numpy.array( ( [266.41683, -29.00781, 2000, 0, 0, 0, 0, 0, 0, 3600, 3600, 0], - [50.67375, -37.20833, 528, 0, 0, 0, 178e6, -0.51, 0, 0, 0, 0], # For - [201.36667, -43.01917, 1370, 0, 0, 0, 200e6, -0.50, 0, 0, 0, 0], # Cen - [139.52500, -12.09556, 280, 0, 0, 0, 200e6, -0.96, 0, 0, 0, 0], # Hyd - [79.95833, -45.77889, 390, 0, 0, 0, 200e6, -0.99, 0, 0, 0, 0], # Pic - [252.78333, 4.99250, 377, 0, 0, 0, 200e6, -1.07, 0, 0, 0, 0], # Her - [187.70417, 12.39111, 861, 0, 0, 0, 200e6, -0.86, 0, 0, 0, 0], # Vir - [83.63333, 22.01444, 1340, 0, 0, 0, 200e6, -0.22, 0, 0, 0, 0], # Tau - [299.86667, 40.73389, 7920, 0, 0, 0, 200e6, -0.78, 0, 0, 0, 0], # Cyg - [350.86667, 58.81167, 11900, 0, 0, 0, 200e6, -0.41, 0, 0, 0, 0], # Cas + [ + 50.67375, + -37.20833, + 528, + 0, + 0, + 0, + 178e6, + -0.51, + 0, + 0, + 0, + 0, + ], # For + [ + 201.36667, + -43.01917, + 1370, + 0, + 0, + 0, + 200e6, + -0.50, + 0, + 0, + 0, + 0, + ], # Cen + [ + 139.52500, + -12.09556, + 280, + 0, + 0, + 0, + 200e6, + -0.96, + 0, + 0, + 0, + 0, + ], # Hyd + [ + 79.95833, + -45.77889, + 390, + 0, + 0, + 0, + 200e6, + -0.99, + 0, + 0, + 0, + 0, + ], # Pic + [ + 252.78333, + 4.99250, + 377, + 0, + 0, + 0, + 200e6, + -1.07, + 0, + 0, + 0, + 0, + ], # Her + [ + 187.70417, + 12.39111, + 861, + 0, + 0, + 0, + 200e6, + -0.86, + 0, + 0, + 0, + 0, + ], # Vir + [ + 83.63333, + 22.01444, + 1340, + 0, + 0, + 0, + 200e6, + -0.22, + 0, + 0, + 0, + 0, + ], # Tau + [ + 299.86667, + 40.73389, + 7920, + 0, + 0, + 0, + 200e6, + -0.78, + 0, + 0, + 0, + 0, + ], # Cyg + [ + 350.86667, + 58.81167, + 11900, + 0, + 0, + 0, + 200e6, + -0.41, + 0, + 0, + 0, + 0, + ], # Cas ) ) @@ -81,11 +198,15 @@ def make_sky_model(sky0, settings, radius_deg, flux_min_outer_jy): sky_outer.num_sources, ) sky_outer.append(sky_inner) - LOG.info("Number of sources in output sky model: %d", sky_outer.num_sources) + LOG.info( + "Number of sources in output sky model: %d", sky_outer.num_sources + ) return sky_outer -def make_plot(prefix, field_name, metric_key, results, axis_length, axis_bandwidth): +def make_plot( + prefix, field_name, metric_key, results, axis_length, axis_bandwidth +): """Plot selected results.""" # Axis setup. ax1 = plt.subplot(111) @@ -143,13 +264,19 @@ def create_imagers_multi_channel( channel_end = channel_start + channel_step - 1 if channel_end >= num_channels - 1: channel_end = num_channels - 1 - key = "length_%.0f_channels_%d-%d" % (length_sec, channel_start, channel_end) + key = "length_%.0f_channels_%d-%d" % ( + length_sec, + channel_start, + channel_end, + ) if key in results: continue imager = oskar.Imager(settings=settings_img) freq_min = freq_start + freq_inc * (channel_start - 0.5) freq_max = freq_start + freq_inc * (channel_end + 0.5) - LOG.info('"%s" has frequency range %.5e to %.5e Hz', key, freq_min, freq_max) + LOG.info( + '"%s" has frequency range %.5e to %.5e Hz', key, freq_min, freq_max + ) imager.channel_snapshots = True imager.grid_on_gpu = True imager.freq_min_hz = freq_min @@ -158,7 +285,9 @@ def create_imagers_multi_channel( imager.time_max_utc = time_max imager.output_root = prefix_field + "_" + key imagers.append((key, imager)) - LOG.info("Will make %d images of %d channels each.", len(imagers), channel_step) + LOG.info( + "Will make %d images of %d channels each.", len(imagers), channel_step + ) return imagers @@ -198,7 +327,9 @@ def read_coords(out0_name, imagers): # Read both visibility blocks. LOG.info("Reading blocks %d/%d", i_block + 1, hdr1.num_blocks) tasks_read.append( - executor.submit(block_diff[i_block % 2].read, hdr1, handle1, i_block) + executor.submit( + block_diff[i_block % 2].read, hdr1, handle1, i_block + ) ) # Ensure reads have finished. @@ -269,7 +400,13 @@ def read_data(out0_name, out_name, prefix_field, imagers): def run_set_new( - base_settings_sim, base_settings_img, sky0, prefix_field, field, length_sec, results + base_settings_sim, + base_settings_img, + sky0, + prefix_field, + field, + length_sec, + results, ): """Runs a set of simulations for one field.""" # Load telescope model. @@ -390,7 +527,10 @@ def run_set( freq_min = freq_mid - (bandwidth * 1e6 + freq_inc) * 0.5 freq_max = freq_mid + (bandwidth * 1e6 + freq_inc) * 0.5 LOG.info( - '"%s" has frequency range %.5e to %.5e Hz', key, freq_min, freq_max + '"%s" has frequency range %.5e to %.5e Hz', + key, + freq_min, + freq_max, ) imager.freq_min_hz = freq_min imager.freq_max_hz = freq_max @@ -417,8 +557,12 @@ def run_set( block1.read(hdr1, handle1, i_block) tasks = [] for i, (key, imager) in enumerate(imagers): - LOG.info('Updating imager %d/%d "%s"...', i + 1, len(imagers), key) - tasks.append(executor.submit(imager.update_from_block, hdr1, block1)) + LOG.info( + 'Updating imager %d/%d "%s"...', i + 1, len(imagers), key + ) + tasks.append( + executor.submit(imager.update_from_block, hdr1, block1) + ) concurrent.futures.wait(tasks) for i, (key, imager) in enumerate(imagers): imager.coords_only = False @@ -440,7 +584,9 @@ def run_set( block1.cross_correlations()[...] -= block2.cross_correlations() tasks = [] for _, imager in imagers: - tasks.append(executor.submit(imager.update_from_block, hdr1, block1)) + tasks.append( + executor.submit(imager.update_from_block, hdr1, block1) + ) concurrent.futures.wait(tasks) executor.shutdown() del handle1, handle2, hdr1, hdr2, block1, block2 @@ -462,10 +608,10 @@ def run_set( "image_medianabs": numpy.median(numpy.abs(image)), "image_mean": numpy.mean(image), "image_std": numpy.std(image), - "image_rms": numpy.sqrt(numpy.mean(image ** 2)), + "image_rms": numpy.sqrt(numpy.mean(image**2)), "image_centre_mean": numpy.mean(centre), "image_centre_std": numpy.std(centre), - "image_centre_rms": numpy.sqrt(numpy.mean(centre ** 2)), + "image_centre_rms": numpy.sqrt(numpy.mean(centre**2)), } del image, output @@ -505,7 +651,9 @@ def main(): } base_settings_img = { "image": { - "double_precision": base_settings_sim["simulator"]["double_precision"], + "double_precision": base_settings_sim["simulator"][ + "double_precision" + ], "size": 8192, "fov_deg": 5.0, "fft/use_gpu": True, diff --git a/tests/rfi/pycraf_scripts/test_pycraf_propagation.py b/tests/rfi/pycraf_scripts/test_pycraf_propagation.py index e4c7b68246e49a41c3a0dfbb5dc3f1a56ba58277..939f27f240dd658c950272b1709513a5d1d0f4e1 100644 --- a/tests/rfi/pycraf_scripts/test_pycraf_propagation.py +++ b/tests/rfi/pycraf_scripts/test_pycraf_propagation.py @@ -2,13 +2,14 @@ Unit tests for pycraf propagation calculation """ -import numpy as np -import pytest import re import unittest -from pycraf import pathprof from urllib.error import HTTPError +import numpy as np +import pytest +from pycraf import pathprof + from rfi.pycraf_scripts.SKA_low_pycraf_propagation import calc_prop_atten @@ -28,7 +29,10 @@ class TestPYprop(unittest.TestCase): temperature = 290.0 # u.K pressure = 1013.0 # u.hPa timepercent = 0.02 # u.percent # see P.452 for explanation - h_tg, h_rg = 175, 2.1 # u.m # height of the receiver and transmitter above gnd + h_tg, h_rg = ( + 175, + 2.1, + ) # u.m # height of the receiver and transmitter above gnd # clutter type for transmitter/receiver zone_t, zone_r = pathprof.CLUTTER.UNKNOWN, pathprof.CLUTTER.UNKNOWN @@ -59,7 +63,7 @@ class TestPYprop(unittest.TestCase): hprof_step = 10000 # u.m # resolution of solution # DTV Bickley site Perth 180MHz Seven - tx_name, tx_lon, tx_lat = "TVW6_26624", 116.061666667, -32.0127777778 + tx_lon, tx_lat = 116.061666667, -32.0127777778 Atten_ant, freqs, apparent_coords, ska_low = calc_prop_atten( freq, @@ -88,7 +92,8 @@ class TestPYprop(unittest.TestCase): # After spending somet tiem on this, I saw that the rows are shifted # in the output array, but I don't know why # original values - # output = np.array([[172.82915997, 172.71465099, 172.56585665], [172.84376699, 172.73321301, 172.5883582], + # output = np.array([[172.82915997, 172.71465099, 172.56585665], + # [172.84376699, 172.73321301, 172.5883582], # [172.83012078, 172.71681939, 172.56922449]]) # for now, I'll update them to pass test output_atten = np.array( diff --git a/tests/rfi/pycraf_scripts/test_rfi_propagation.py b/tests/rfi/pycraf_scripts/test_rfi_propagation.py index 5a01f23f718a06a3da469fb95fc4f6924b41f6e8..26e5b5922f9424c5e6df2423e5e422a86cc590c1 100644 --- a/tests/rfi/pycraf_scripts/test_rfi_propagation.py +++ b/tests/rfi/pycraf_scripts/test_rfi_propagation.py @@ -1,19 +1,18 @@ +import os import re +from unittest.mock import Mock, patch from urllib.error import HTTPError import numpy as np -import os -from unittest.mock import patch, Mock - import pytest from rfi.pycraf_scripts.SKA_low_RFI_propagation import ( - rfi_attenuation, + calculate_apparent_power, cli_parser, - transmitter_below_horizon, construct_csv_file_path, main, - calculate_apparent_power, + rfi_attenuation, + transmitter_below_horizon, ) TEXT_TO_REPLACE_IN_PATH = "tests/rfi/pycraf_scripts/test_rfi_propagation.py" @@ -26,9 +25,11 @@ OUTPUTDIR = "/where/file/should/go/" ) def test_transmitter_below_horizon(check_horizon, elevation, expected_result): """ - If check_horizon is "True", when the elevation of a transmitter is below zero, - the function returns False, else returns True. - If check_horizon is "False", function returns False regardless of transmitter elevation. + If check_horizon is "True", when the elevation of + a transmitter is below zero, the function returns False, + else returns True. + If check_horizon is "False", function returns False + regardless of transmitter elevation. Note: check_horizon has to be a string. """ @@ -39,16 +40,23 @@ def test_transmitter_below_horizon(check_horizon, elevation, expected_result): @pytest.mark.parametrize( "file_name, filtered_transmitter, expected_csv_name", [ - # if file_name is "infile", generate file name from existing filtered_transmitter string - ("infile", "transmitter_file_name", f"{OUTPUTDIR}transmitter_file__el.csv"), + # if file_name is "infile", generate file name + # from existing filtered_transmitter string + ( + "infile", + "transmitter_file_name", + f"{OUTPUTDIR}transmitter_file__el.csv", + ), ("infile", "dir/tr_files", f"{OUTPUTDIR}tr_f_el.csv"), - # if file_name also contains path, append that to full path of working directory + # if file_name also contains path, append that to + # full path of working directory ( "dir/my-transm-file.csv", "existing_transmitter_file", os.getcwd() + "/dir/my-transm-file.csv", ), - # if file_name is not a path, the resulting path will include the outputdir value + # if file_name is not a path, the resulting path + # will include the outputdir value ( "want_my_file_to_be_called_this", "some_existing_file_with_transmitters", @@ -56,8 +64,12 @@ def test_transmitter_below_horizon(check_horizon, elevation, expected_result): ), ], ) -def test_construct_csv_file_path(file_name, filtered_transmitter, expected_csv_name): - result = construct_csv_file_path(file_name, filtered_transmitter, OUTPUTDIR) +def test_construct_csv_file_path( + file_name, filtered_transmitter, expected_csv_name +): + result = construct_csv_file_path( + file_name, filtered_transmitter, OUTPUTDIR + ) assert result == expected_csv_name @@ -70,10 +82,13 @@ def test_calculate_apparent_power(): transm_power = 4000 # W attenuat = np.ones((3, 4)) # 3 antennas, 4 channels - result = calculate_apparent_power(ntimes, frequency_range, transm_power, attenuat) + result = calculate_apparent_power( + ntimes, frequency_range, transm_power, attenuat + ) assert result.shape == (2, 3, 4) # ntimes x antennas x channels - # because attenuat is constant at all points, this value is the same everywhere in result + # because attenuat is constant at all points, + # this value is the same everywhere in result assert (result.round(5) == 0.01783 + 0j).all() @@ -95,7 +110,8 @@ def test_calculate_apparent_power_time_var(): assert result.shape == (2, 3, 4) # ntimes x antennas x channels # first time sample, same for every frequency and station assert (abs(result[0]).round(5) == 0.00737).all() - # second time sample, diff from first, but same for every frequency and station + # second time sample, diff from first, but same for + # every frequency and station assert (abs(result[1]).round(5) == 0.01598).all() @@ -111,7 +127,11 @@ def test_calculate_apparent_power_freq_var(): attenuat = np.ones((3, 4)) # 3 antennas, 4 channels result = calculate_apparent_power( - ntimes, frequency_range, transm_power, attenuat, frequency_variable=True + ntimes, + frequency_range, + transm_power, + attenuat, + frequency_variable=True, ) assert result.shape == (2, 3, 4) # ntimes x antennas x channels @@ -148,8 +168,9 @@ def test_calculate_apparent_power_time_freq_var(): assert result.shape == (2, 3, 4) # ntimes x antennas x channels # 0th station (for any station the results are the same) - # all of the values are time and frequency variable, hence they are all unique - # the length of unique values in the array matches the length of the flattened array + # all of the values are time and frequency variable, + # hence they are all unique the length of unique values in + # the array matches the length of the flattened array assert len(np.unique(result[:, 0, :])) == len(result[:, 0, :].flatten()) @@ -159,13 +180,16 @@ def test_calculate_apparent_power_time_freq_var(): @patch("rfi.pycraf_scripts.SKA_low_RFI_propagation.write_settings", Mock()) def test_rfi_attenuation(mock_to_csv): """ - With the given frequency range, and the transmitter central frequency from the - transmitter file, the simulation will produce emitted powers for two frequency channels, + With the given frequency range, and the transmitter + central frequency from the transmitter file, the simulation + will produce emitted powers for two frequency channels, which we test here. """ file_path = __file__ workbook_location = re.sub( - TEXT_TO_REPLACE_IN_PATH, "rfi/test_data/SKA_Low_test_REV3.txt", file_path + TEXT_TO_REPLACE_IN_PATH, + "rfi/test_data/SKA_Low_test_REV3.txt", + file_path, ) srtm_dir = re.sub(TEXT_TO_REPLACE_IN_PATH, "rfi/test_data", file_path) transmitter_file = re.sub( @@ -210,7 +234,9 @@ def test_rfi_attenuation(mock_to_csv): assert round(result_power[0, 1, 1], 3) == round(expected_power[0, 1, 1], 3) result_freq = result[0].rfi_signal.frequency # in Hz - assert (np.round(result_freq, 3) == np.round(np.array([1.75e08, 1.80e08]), 3)).all() + assert ( + np.round(result_freq, 3) == np.round(np.array([1.75e08, 1.80e08]), 3) + ).all() @pytest.mark.xfail(raises=HTTPError, reason="SMRT data website is down") @@ -220,17 +246,22 @@ def test_rfi_attenuation(mock_to_csv): def test_rfi_attenuation_bug_fix(mock_to_csv): """ Bug Fix: - When multiple transmitters are used, and the first one is skipped due to its central frequency - not fitting within the simulated frequency range, the second transmitter is throwing a TypeError - because start and end frequencies of the range are not reset to the original value but left as None. - - We have 2 transmitters. One is out of frequency range. The second is usable and calculations - should be performed for it, inc. saving outputs. + When multiple transmitters are used, and the first one is + skipped due to its central frequency not fitting within + the simulated frequency range, the second transmitter is + throwing a TypeError because start and end frequencies of + the range are not reset to the original value but left as None. + + We have 2 transmitters. One is out of frequency range. + The second is usable and calculations should be + performed for it, inc. saving outputs. """ file_path = __file__ workbook_location = re.sub( - TEXT_TO_REPLACE_IN_PATH, "rfi/test_data/SKA_Low_test_REV3.txt", file_path + TEXT_TO_REPLACE_IN_PATH, + "rfi/test_data/SKA_Low_test_REV3.txt", + file_path, ) srtm_dir = re.sub(TEXT_TO_REPLACE_IN_PATH, "rfi/test_data", file_path) transmitter_file = re.sub( @@ -267,7 +298,10 @@ def test_rfi_attenuation_bug_fix(mock_to_csv): @patch("rfi.pycraf_scripts.SKA_low_RFI_propagation.check_for_old_run", Mock()) @patch("rfi.pycraf_scripts.SKA_low_RFI_propagation.write_settings", Mock()) -@patch("rfi.pycraf_scripts.SKA_low_RFI_propagation.DataCube.export_to_hdf5", Mock()) +@patch( + "rfi.pycraf_scripts.SKA_low_RFI_propagation.DataCube.export_to_hdf5", + Mock(), +) @patch("rfi.pycraf_scripts.SKA_low_RFI_propagation.cli_parser") @patch("rfi.pycraf_scripts.SKA_low_RFI_propagation.rfi_attenuation") def test_rfi_propagation_main(new_rfi_atten, mock_parser): @@ -276,7 +310,9 @@ def test_rfi_propagation_main(new_rfi_atten, mock_parser): """ file_path = __file__ workbook_location = re.sub( - TEXT_TO_REPLACE_IN_PATH, "rfi/test_data/SKA_Low_test_REV3.txt", file_path + TEXT_TO_REPLACE_IN_PATH, + "rfi/test_data/SKA_Low_test_REV3.txt", + file_path, ) srtm_dir = re.sub(TEXT_TO_REPLACE_IN_PATH, "rfi/test_data", file_path) transmitter_file = re.sub( @@ -304,7 +340,8 @@ def test_rfi_propagation_main(new_rfi_atten, mock_parser): mock_parser().parse_args.return_value = args - # replace the run within main with this. This way the optional arguments can be passed in. + # replace the run within main with this. + # This way the optional arguments can be passed in. new_rfi_atten.return_value = rfi_attenuation( args, save_azel_to_file=False, hdf5_only=True ) diff --git a/tests/rfi/rascil_scripts/test_regression_rascil_simulate_low_rfi_visibility_propagation.py b/tests/rfi/rascil_scripts/test_regression_rascil_simulate_low_rfi_visibility_propagation.py index d2ce337deebe777f25c29b761abdf3de71725e6d..d1be8236c2bdc35a7da27da27cba521a2f2d2df8 100644 --- a/tests/rfi/rascil_scripts/test_regression_rascil_simulate_low_rfi_visibility_propagation.py +++ b/tests/rfi/rascil_scripts/test_regression_rascil_simulate_low_rfi_visibility_propagation.py @@ -6,21 +6,23 @@ import logging import os import re import tempfile +from unittest.mock import Mock, patch import numpy import pytest -from unittest.mock import patch, Mock - from rascil.apps.rascil_imager import cli_parser as imager_cli_parser from rascil.apps.rascil_imager import imager from rascil.apps.rascil_vis_ms import cli_parser as vis_cli_parser from rascil.apps.rascil_vis_ms import visualise -from rfi.rascil_scripts.power_spectrum import cli_parser as powerspectrum_cli_parser -from rfi.rascil_scripts.power_spectrum import power_spectrum from rascil.processing_components import import_image_from_fits, qa_image + +from rfi.rascil_scripts.power_spectrum import ( + cli_parser as powerspectrum_cli_parser, +) +from rfi.rascil_scripts.power_spectrum import power_spectrum from rfi.rascil_scripts.simulate_low_rfi_visibility_propagation import ( - rfi_simulation, cli_parser, + rfi_simulation, ) log = logging.getLogger("rascil-logger") @@ -52,28 +54,34 @@ PERSIST = os.getenv("RASCIL_PERSIST", False) ], ) @patch( - "rfi.rascil_scripts.simulate_low_rfi_visibility_propagation.check_for_old_run", + "rfi.rascil_scripts.simulate_low_rfi_visibility" + "_propagation.check_for_old_run", Mock(), ) @patch( - "rfi.rascil_scripts.simulate_low_rfi_visibility_propagation.write_settings", Mock() + "rfi.rascil_scripts.simulate_low_rfi_visibility" + "_propagation.write_settings", + Mock(), ) def test_rfi_simulation(use_dask, msout): file_path = __file__ input_file = re.sub( - "tests/rfi/rascil_scripts/test_regression_rascil_simulate_low_rfi_visibility_propagation.py", + "tests/rfi/rascil_scripts/test_regression_rascil_" + "simulate_low_rfi_visibility_propagation.py", "rfi/test_data/tv_transmitter_attenuation_cube.hdf5", file_path, ) antenna_file = re.sub( - "tests/rfi/rascil_scripts/test_regression_rascil_simulate_low_rfi_visibility_propagation.py", - "rfi/data/telescope_files/SKA1-LOW_SKO-0000422_Rev3_38m_SKALA4_spot_frequencies.tm/layout_wgs84.txt", + "tests/rfi/rascil_scripts/test_regression_rascil_" + "simulate_low_rfi_visibility_propagation.py", + "rfi/data/telescope_files/SKA1-LOW_SKO-0000422_Rev3_" + "38m_SKALA4_spot_frequencies.tm/layout_wgs84.txt", file_path, ) - expected_image_max = 20438743.228243344 - expected_image_min = -10312773.328625064 - expected_image_rms = 127784.79100904595 + expected_image_max = 20438501.83649988 + expected_image_min = -10312397.909026483 + expected_image_rms = 127783.89983864271 with tempfile.TemporaryDirectory() as tempdirname: @@ -144,7 +152,12 @@ def test_rfi_simulation(use_dask, msout): dirty = import_image_from_fits(dirtyname) qa = qa_image(dirty) - assert dirty["pixels"].shape == (1, 1, args.imaging_npixel, args.imaging_npixel) + assert dirty["pixels"].shape == ( + 1, + 1, + args.imaging_npixel, + args.imaging_npixel, + ) numpy.testing.assert_allclose( qa.data["max"], expected_image_max, atol=1e-7, rtol=1e-7 ) diff --git a/tests/rfi/rascil_scripts/test_unit_simulate_low_rfi_visibility_prop.py b/tests/rfi/rascil_scripts/test_unit_simulate_low_rfi_visibility_prop.py index ab51ba1e828f663526c409c326898778ea6ad177..03a9c8c70a827fb8d415c38deb98d3358daced5e 100644 --- a/tests/rfi/rascil_scripts/test_unit_simulate_low_rfi_visibility_prop.py +++ b/tests/rfi/rascil_scripts/test_unit_simulate_low_rfi_visibility_prop.py @@ -2,14 +2,13 @@ import os import tempfile from unittest.mock import Mock, patch -import numpy as np import astropy.units as u - +import numpy as np from astropy.coordinates import SkyCoord from rascil.data_models import PolarisationFrame from rascil.processing_components import ( - create_named_configuration, create_blockvisibility, + create_named_configuration, ) from rfi.rascil_scripts.simulate_low_rfi_visibility_propagation import ( @@ -72,7 +71,9 @@ def test_load_beam_gain_and_ra_dec_dont_use_beam(): assert result[2] == -45.0 # DEC -@patch("rfi.rascil_scripts.simulate_low_rfi_visibility_propagation.numpy.loadtxt") +@patch( + "rfi.rascil_scripts.simulate_low_rfi_visibility_propagation.numpy.loadtxt" +) def test_load_beam_gain_and_ra_dec_numpy_beam(mock_load): """ args.use_beamgain = True @@ -95,7 +96,9 @@ def test_load_beam_gain_and_ra_dec_numpy_beam(mock_load): ] # 2 sources, 3 frequency channels expected_gain = np.array([[[9, 4, 16]], [[4, 25, 9]]]) - result = load_beam_gain_and_ra_dec(args, source_ids, frequency_channels, Mock()) + result = load_beam_gain_and_ra_dec( + args, source_ids, frequency_channels, Mock() + ) assert (result[0] == expected_gain).all() # beam gain assert result[1] == 0.0 # RA @@ -117,7 +120,8 @@ def test_load_beam_gain_and_ra_dec_hdf_beam(): test_beam_gain_cube = BeamGainDataCube( ra, dec, - "01-01-2000 09:32:00.000", # time of observation, same as one in OSKAR settings + # time of observation, same as one in OSKAR settings + "01-01-2000 09:32:00.000", frequency_channels, source_ids, nstations, # nstations @@ -148,7 +152,9 @@ def test_load_beam_gain_and_ra_dec_hdf_beam(): ) test_beam_gain_cube.export_to_hdf5(hdf_file) - result = load_beam_gain_and_ra_dec(args, source_ids, frequency_channels, Mock()) + result = load_beam_gain_and_ra_dec( + args, source_ids, frequency_channels, Mock() + ) os.remove(hdf_file) @@ -169,10 +175,13 @@ def test_get_chunk_start_times(): nintegr_per_chunk = 6 int_time = 10.0 # [s] - result = get_chunk_start_times(time_range, int_time, nintegr_per_chunk, Mock()) + result = get_chunk_start_times( + time_range, int_time, nintegr_per_chunk, Mock() + ) assert len(result) == 60 - assert result[0] == np.array([0]) # start time of first chunk, 0 seconds - assert result[-1] == np.array( - [3540] - ) # start time of last chunk, 3540 seconds, 60 s before end of time_range (1h) + # start time of first chunk, 0 seconds + assert result[0] == np.array([0]) + # start time of last chunk, 3540 seconds, + # 60 s before end of time_range (1h) + assert result[-1] == np.array([3540]) diff --git a/tests/rfi/rfi_interface/conftest.py b/tests/rfi/rfi_interface/conftest.py index 6ae7b03bfcb4e0fdb8e560628a79eaff49cc3b95..4263d62ea88e2c021613b49f1682494aa3c9c198 100644 --- a/tests/rfi/rfi_interface/conftest.py +++ b/tests/rfi/rfi_interface/conftest.py @@ -2,9 +2,9 @@ import numpy as np import pytest from rfi.rfi_interface.rfi_data_cube import ( - RFISource, - RFISignal, DataCubePerSource, + RFISignal, + RFISource, ) @@ -31,11 +31,19 @@ def rfi_signal(): "station_id": np.array(["ska_low_station1", "ska_low_station2"]), "source_id": "test", "time": np.array( - ["2021-04-11 01:20:55", "2021-04-11 01:30:55", "2021-04-11 01:40:55"] + [ + "2021-04-11 01:20:55", + "2021-04-11 01:30:55", + "2021-04-11 01:40:55", + ] ), "frequency": np.array([250.1 * 1e6, 250.2 * 1e6, 250.3 * 1e6]), - "azimuth": np.array([[12.0, 11.0] * 3]), # 3 for the three time-samples - "elevation": np.array([[0.5, 0.6] * 3]), # 3 for the three time-samples + "azimuth": np.array( + [[12.0, 11.0] * 3] + ), # 3 for the three time-samples + "elevation": np.array( + [[0.5, 0.6] * 3] + ), # 3 for the three time-samples "distance": np.array( [[150.0 * 1e3, 160.0 * 1e3] * 3] ), # 3 for the three time-samples diff --git a/tests/rfi/rfi_interface/test_rfi_data_cube.py b/tests/rfi/rfi_interface/test_rfi_data_cube.py index 79c448143b37fd21f4a2829d3262c76a9183694f..5f5da753922a9320044ec8b870988c7521181403 100644 --- a/tests/rfi/rfi_interface/test_rfi_data_cube.py +++ b/tests/rfi/rfi_interface/test_rfi_data_cube.py @@ -1,15 +1,14 @@ from copy import deepcopy -from unittest.mock import patch, Mock +from dataclasses import fields +from unittest.mock import Mock, patch -import pytest import numpy as np - -from dataclasses import fields +import pytest from rfi.rfi_interface.rfi_data_cube import ( + BeamGainDataCube, DataCube, DataCubePerSource, - BeamGainDataCube, ) @@ -33,7 +32,10 @@ def test_rfi_signal(rfi_signal): """ assert ( np.column_stack( - (rfi_signal.frequency, np.column_stack(rfi_signal.apparent_power[0, :])) + ( + rfi_signal.frequency, + np.column_stack(rfi_signal.apparent_power[0, :]), + ) ) == np.array( [ @@ -47,13 +49,15 @@ def test_rfi_signal(rfi_signal): rfi_signal.azimuth == np.array([[12.0, 11.0] * 3]) ).all() # there are three time-samples assert ( - rfi_signal.station_id == np.array(["ska_low_station1", "ska_low_station2"]) + rfi_signal.station_id + == np.array(["ska_low_station1", "ska_low_station2"]) ).all() def test_data_cube_per_source(data_cube_per_source, rfi_signal): """ - DataCubePerSource object is created with the correct values from input dict. + DataCubePerSource object is created with the + correct values from input dict. :param data_cube_per_source: DataCubePerSource fixture from conftest.py :param rfi_signal: RFISignal fixture from conftest.py @@ -63,12 +67,15 @@ def test_data_cube_per_source(data_cube_per_source, rfi_signal): assert data_cube_per_source.rfi_signal == rfi_signal -def test_data_cube_per_source_wrong_source_id_in_signal(rfi_signal, rfi_source): +def test_data_cube_per_source_wrong_source_id_in_signal( + rfi_signal, rfi_source +): """ The id of the RFISource object doesn't match the source id in any of the RFISignal objects within the DataCubePerSource object. - The source id in rfi_source, source_at_ska and rfi_signal is "test" by default. + The source id in rfi_source, source_at_ska and rfi_signal + is "test" by default. :param source_at_ska: RFISourceAtSKA fixture from conftest.py :param rfi_signal: RFISignal fixture from conftest.py @@ -142,7 +149,8 @@ def test_data_cube_append(data_cube_per_source): @pytest.mark.parametrize("source_id", [None, "", b""]) def test_data_cube_error_missing_source_id(source_id, data_cube_per_source): """ - Raise ValueError if the data being appended to DataCube does not contain a valid source_id. + Raise ValueError if the data being appended to + DataCube does not contain a valid source_id. :param data_cube_per_source: fixtures from conftest.py """ @@ -161,7 +169,8 @@ def test_data_cube_error_missing_source_id(source_id, data_cube_per_source): def test_data_cube_error_wrong_times(data_cube_per_source): """ - Data appended to DataCube has to refer to the same time samples as the DataCube was + Data appended to DataCube has to refer to the + same time samples as the DataCube was initialized with; otherwise raise ValueError :param data_cube_per_source: fixtures from conftest.py @@ -180,7 +189,8 @@ def test_data_cube_error_wrong_times(data_cube_per_source): def test_data_cube_error_wrong_frequency(data_cube_per_source): """ - Data appended to DataCube has to refer to the same frequency samples as the DataCube was + Data appended to DataCube has to refer to the + same frequency samples as the DataCube was initialized with; otherwise raise ValueError :param data_cube_per_source: fixtures from conftest.py @@ -199,7 +209,8 @@ def test_data_cube_error_wrong_frequency(data_cube_per_source): def test_data_cube_error_wrong_station_id(data_cube_per_source): """ - Data appended to DataCube has to refer to the same station ids as the DataCube was + Data appended to DataCube has to refer to the + same station ids as the DataCube was initialized with; otherwise raise ValueError :param data_cube_per_source: fixtures from conftest.py @@ -229,7 +240,9 @@ def beam_gain_cube(): rfi_source_ids = np.array([b"S1", b"S2"]) nstations = 5 - result = BeamGainDataCube(ra, dec, obs_time, freq_chans, rfi_source_ids, nstations) + result = BeamGainDataCube( + ra, dec, obs_time, freq_chans, rfi_source_ids, nstations + ) yield result diff --git a/tests/rfi/rfi_interface/test_rfi_source_signal_interface.py b/tests/rfi/rfi_interface/test_rfi_source_signal_interface.py index 78e70bb359936747525860b8765bdbaf3e902378..554ea91941e1da9c62f59ab4f5628bd9679ae0aa 100644 --- a/tests/rfi/rfi_interface/test_rfi_source_signal_interface.py +++ b/tests/rfi/rfi_interface/test_rfi_source_signal_interface.py @@ -1,14 +1,14 @@ -import pytest +from unittest.mock import MagicMock, patch -from unittest.mock import patch, MagicMock +import pytest from rfi.rfi_interface.rfi_source_signal_interface import ( - rfi_interface, - FREQ_RANGE, - N_CHANNELS, + FREQ_END, FREQ_KEYS, + FREQ_RANGE, FREQ_START, - FREQ_END, + N_CHANNELS, + rfi_interface, validate_frequency, ) @@ -54,7 +54,8 @@ def test_initialize_data_cube_error_two_given(): def test_initialize_data_cube(): """ - If all of FREQ_KEYS are provided, add their formatted version to the input_dictionary + If all of FREQ_KEYS are provided, add their + formatted version to the input_dictionary """ input_args = { f"<{N_CHANNELS}>": 5, @@ -98,7 +99,10 @@ def test_rfi_interface_tv_antenna(mock_cli_parser, mock_rfi_atten): Happy path for rfi_interface, with source_type = tv_antenna """ source_type = "tv_antenna" - input_args = {"--frequency_range": ["170.5", "184.5"], "--n_time_chunks": "1"} + input_args = { + "--frequency_range": ["170.5", "184.5"], + "--n_time_chunks": "1", + } mock_cli_parser.return_value = MagicMock(n_time_chunks="1") mock_rfi_atten.return_value = MagicMock() diff --git a/tests/rfi/test_common_tools.py b/tests/rfi/test_common_tools.py index 52886d40815465f11125e4cdc1adc72a0dae3939..0adac04c7c7eaead0ba71d17a739f8264871363c 100644 --- a/tests/rfi/test_common_tools.py +++ b/tests/rfi/test_common_tools.py @@ -2,7 +2,10 @@ import numpy as np import pandas as pd import pytest -from rfi.common_tools import calculate_bandwidth_channel_range, default_or_antenna_value +from rfi.common_tools import ( + calculate_bandwidth_channel_range, + default_or_antenna_value, +) def test_calculate_bandwidth_channel_range(): @@ -19,7 +22,9 @@ def test_calculate_bandwidth_channel_range(): tx_central_freq = 150 freq_array = np.linspace(50, 250, 10) - result = calculate_bandwidth_channel_range(bandwidth, freq_array, tx_central_freq) + result = calculate_bandwidth_channel_range( + bandwidth, freq_array, tx_central_freq + ) assert result[0] == 2 # number of channels assert round(result[1], 2) == 138.89 # start frequency @@ -28,7 +33,8 @@ def test_calculate_bandwidth_channel_range(): def test_calculate_bandwidth_channel_range_index_error(): """ - IndexError caused by dtv_range calculated from bandwidth and central_frequency + IndexError caused by dtv_range calculated + from bandwidth and central_frequency not matching the freq_array range. """ bandwidth = 7.0 @@ -39,7 +45,9 @@ def test_calculate_bandwidth_channel_range_index_error(): freq_end = 210.0 freq_array = np.linspace(freq_start, freq_end, n_freqs) - result = calculate_bandwidth_channel_range(bandwidth, freq_array, tx_central_freq) + result = calculate_bandwidth_channel_range( + bandwidth, freq_array, tx_central_freq + ) assert result == (None, None, None, None) @@ -48,7 +56,12 @@ def test_calculate_bandwidth_channel_range_index_error(): "property_name, set_value, default_value, expected_value", [ ("central frequency", "True", 10.0, 10.0), # return default value - ("central frequency", "False", 10.0, 180.0), # return value from antenna_data + ( + "central frequency", + "False", + 10.0, + 180.0, + ), # return value from antenna_data ( "bandwidth", "False", @@ -79,10 +92,13 @@ def test_default_or_antenna_value_wrong_property(): antenna_data = pd.Series() with pytest.raises(KeyError) as e: - default_or_antenna_value(property_name, set_value, default_value, antenna_data) + default_or_antenna_value( + property_name, set_value, default_value, antenna_data + ) assert ( - e.value.args[0] == "The property you provided (some-property) is not valid. " + e.value.args[0] + == "The property you provided (some-property) is not valid. " "Currently only the following can be used: " "dict_keys(['bandwidth', 'central frequency'])" ) diff --git a/tests/rfi/test_old_run_checker.py b/tests/rfi/test_old_run_checker.py index 81fa49df1ae36d331d88c0f42c0ac952ab22df38..1b004c3c99fea91e9f9f5c1765e2a8b93257442d 100644 --- a/tests/rfi/test_old_run_checker.py +++ b/tests/rfi/test_old_run_checker.py @@ -1,10 +1,10 @@ +import csv import logging -import unittest import os -import csv import tracemalloc +import unittest +from unittest.mock import call, patch -from unittest.mock import patch, call from rfi.old_run_checker import check_for_old_run, write_settings tracemalloc.start() @@ -18,7 +18,12 @@ class TestOldRunChecker(unittest.TestCase): self._settings_file = "test_settings.csv" # Mock settings dict containg basic types. - self._args = {"a": "a/path/to/a/directory", "b": 123.4, "c": False, "d": 2} + self._args = { + "a": "a/path/to/a/directory", + "b": 123.4, + "c": False, + "d": 2, + } self._args_change = { "a": "a/path/to/another/directory", @@ -61,11 +66,16 @@ class TestOldRunChecker(unittest.TestCase): def test_old_run_checker_with_change(self): with patch("logging.Logger.info") as mock_log: check_for_old_run(self._args_change, self._settings_file) - mock_log.assert_has_calls([call("Old settings do not match new settings.")]) + mock_log.assert_has_calls( + [call("Old settings do not match new settings.")] + ) def test_old_run_checker_with_change_order(self): self.assertRaises( - SystemExit, check_for_old_run, self._args_change_order, self._settings_file + SystemExit, + check_for_old_run, + self._args_change_order, + self._settings_file, ) diff --git a/tests/rfi/test_oskar_sim_beam_gain.py b/tests/rfi/test_oskar_sim_beam_gain.py index 970ff895d180887f32152da9e39057418ff67787..4b935014cf661394317d010842adf13c15b36f99 100644 --- a/tests/rfi/test_oskar_sim_beam_gain.py +++ b/tests/rfi/test_oskar_sim_beam_gain.py @@ -1,14 +1,15 @@ import os import re -import pytest +from unittest.mock import Mock, mock_open, patch + import pandas as pd -from unittest.mock import patch, Mock, mock_open +import pytest from rfi.oskar_sim_beam_gain_sing import ( - run_oskar_beam, - cli_parser, _generate_default_oskar_settings, _generate_oskar_settings_from_hdf5, + cli_parser, + run_oskar_beam, run_with_hdf5, ) @@ -39,11 +40,14 @@ ROOT_PATH = "rfi.oskar_sim_beam_gain_sing" @patch(f"{ROOT_PATH}.check_for_old_run", Mock()) @patch(f"{ROOT_PATH}.write_settings", Mock()) @patch(f"{ROOT_PATH}.subprocess.call") -def test_run_oskar_beam(mock_subprocess_call, singularity_exists, expected_app): +def test_run_oskar_beam( + mock_subprocess_call, singularity_exists, expected_app +): """ When oskar_path doesn't exist, subprocess.call is called with the non-singularity version of OSKAR. - When oskar_path exists, subprocess.call is called with the arguments for singularity. + When oskar_path exists, subprocess.call is called + with the arguments for singularity. """ file_path = __file__ transmitter_file = re.sub( @@ -84,7 +88,8 @@ def test_run_oskar_beam(mock_subprocess_call, singularity_exists, expected_app): Mock(side_effect=[True, singularity_exists]), ): - # path the part where we are writing into the oskar_settings.ini file; we don't want that from a test + # path the part where we are writing into the oskar_settings.ini file; + # we don't want that from a test with patch(f"{ROOT_PATH}.open", mock_open()): run_oskar_beam(args) mock_subprocess_call.assert_called_with( @@ -107,7 +112,7 @@ def test_generate_default_oskar_settings( --set_bandwidth = False will take the bandwidth value from the provided dtv_antenna Series. 2. --choose-range = False and --set_bandwidth = True (with the given args) will enter the path where freq_start, freq_end, n_freq, freq_inc are recalculated - """ + """ # noqa: E501 arg_parser = cli_parser() args = arg_parser.parse_args( @@ -127,7 +132,11 @@ def test_generate_default_oskar_settings( ) dtv_antenna = pd.Series( - {"Name": "My fancy antenna", "Frequency(MHz)": 177.5, "Bandwidth(MHz)": 30.0} + { + "Name": "My fancy antenna", + "Frequency(MHz)": 177.5, + "Bandwidth(MHz)": 30.0, + } ) input_dir = "path/to-in-dir" output_dir = "path/to-out-dir" @@ -155,11 +164,12 @@ def test_generate_default_oskar_settings( def test_generate_default_oskar_settings_no_freq_range(): """ - Test that _generate_default_oskar_settings returns None, when the bandwidth-dependent - frequency range cannot be calculated. + Test that _generate_default_oskar_settings returns None, when + the bandwidth-dependent frequency range cannot be calculated. --frequency_range arg does not contain the range provided by - ["Frequency(MHz)" - "Bandwidth(MHz)", "Frequency(MHz)" + "Bandwidth(MHz)"] + ["Frequency(MHz)" - "Bandwidth(MHz)", + "Frequency(MHz)" + "Bandwidth(MHz)"] """ arg_parser = cli_parser() @@ -229,7 +239,9 @@ def test_generate_oskar_settings_hdf5(): assert result["observation"]["num_channels"] == 4 assert result["observation"]["frequency_inc_hz"] == 0.5e6 assert result["beam_pattern"]["sky_model/file"] == "tmp_in.txt" - assert result["beam_pattern"]["station_outputs/text_file/auto_power"] is True + assert ( + result["beam_pattern"]["station_outputs/text_file/auto_power"] is True + ) @patch(f"{ROOT_PATH}.BeamGainDataCube.export_to_hdf5", Mock()) @@ -253,11 +265,17 @@ def test_run_with_hdf5(mock_subprocess_call): These will be automatically removed by code unless test fails. *args needed because the func this one replaces takes an input """ - to_create = ["tmp_out_I.txt", "tmp_out_U.txt", "tmp_out_V.txt", "tmp_out_Q.txt"] + to_create = [ + "tmp_out_I.txt", + "tmp_out_U.txt", + "tmp_out_V.txt", + "tmp_out_Q.txt", + ] for name in to_create: with open(name, "w") as f: f.write("# some comment\n") - # add one line of "beam gain", because there is one channel in the test hdf5 file. + # add one line of "beam gain", because there + # is one channel in the test hdf5 file. f.write("0.0033\n") return 0 @@ -285,7 +303,12 @@ def test_run_with_hdf5(mock_subprocess_call): # code runs without any errors run_with_hdf5( - ["oskar"], args, in_dir, in_dir, telescope_model_name, telescope_model_path + ["oskar"], + args, + in_dir, + in_dir, + telescope_model_name, + telescope_model_path, ) os.remove("oskar_settings.ini")