From 918d42a06a9ab734467808456362d959373913d5 Mon Sep 17 00:00:00 2001 From: mancini Date: Tue, 16 Jun 2020 15:59:23 +0200 Subject: [PATCH] Add self cal algorithm and example script --- lofarimaging/lofarimaging.py | 74 ++++++++++++++++++++++++++-- lofarimaging/singlestationutil.py | 80 ++++++++++++++++++------------- self_cal_demo.py | 66 +++++++++++++++++++++++++ 3 files changed, 184 insertions(+), 36 deletions(-) create mode 100644 self_cal_demo.py diff --git a/lofarimaging/lofarimaging.py b/lofarimaging/lofarimaging.py index 354e619..f5f335d 100644 --- a/lofarimaging/lofarimaging.py +++ b/lofarimaging/lofarimaging.py @@ -1,12 +1,14 @@ """Functions for working with LOFAR single station data""" +import logging from typing import Dict, List -import numpy as np -from numpy.linalg import norm, lstsq -import numexpr as ne + import numba +import numexpr as ne +import numpy as np +import scipy.optimize from astropy.coordinates import SkyCoord, SkyOffsetFrame, CartesianRepresentation - +from numpy.linalg import norm, lstsq __all__ = ["nearfield_imager", "sky_imager", "ground_imager", "skycoord_to_lmn", "calibrate", "simulate_sky_source", "subtract_sources"] @@ -170,6 +172,70 @@ def calibrate(vis, modelvis, maxiter=30, amplitudeonly=True): return residual, gains +def compute_calibrated_model(vis, model_vis, maxiter=30): + n_ant = vis.shape[0] + gain_phase = np.zeros((n_ant), dtype=complex) + + def gains_difference(gain_phase, vis, model_vis): + gains = np.exp(1.j * gain_phase) + gains_mat = np.diag(gains) + gains_mat_star = np.diag(np.conj(gains)) + predicted = gains_mat_star @ model_vis @ gains_mat + residual = np.sum(np.abs(vis - predicted)) + return residual + + result = scipy.optimize.minimize(gains_difference, gain_phase, args=(vis, model_vis), options={'maxiter': maxiter}) + gain_phase = result.x + gains_mat = np.diag(np.exp(1.j * gain_phase)) + + calibrated = np.conj(gains_mat) @ model_vis @ gains_mat + + return calibrated, gains_difference(gain_phase, vis, model_vis) + + +def compute_pointing_matrix(sources_positions, baselines, frequency): + n_baselines = baselines.shape[0] + n_sources = sources_positions.shape[0] + + pointing_matrix = np.zeros(shape=(n_baselines, n_baselines, n_sources), dtype=complex) + + for q in range(n_sources): + pointing_matrix[:, :, q] = simulate_sky_source(sources_positions[q, :], baselines, frequency) + return pointing_matrix + + +def estimate_sources_flux(visibilities, pointing_matrix): + def residual_flux(signal, vis, point_matrix): + residual = vis - point_matrix @ signal + residual = np.sum(np.abs(residual)) + return residual + + n_antennas = visibilities.shape[0] + n_sources = pointing_matrix.shape[2] + lin_visibilities = visibilities.reshape((n_antennas * n_antennas)) + lin_pointing_matrix = pointing_matrix.reshape((n_antennas * n_antennas, n_sources)) + fluxes, *_ = np.linalg.lstsq(lin_pointing_matrix, lin_visibilities) + + result = scipy.optimize.minimize(residual_flux, fluxes, args=(visibilities, pointing_matrix)) + return result.x + + +def estimate_model_visibilities(sources_positions, visibilities, baselines, frequency): + pointing_matrix = compute_pointing_matrix(sources_positions, baselines, frequency) + fluxes = estimate_sources_flux(visibilities, pointing_matrix) + model = pointing_matrix @ fluxes + + return model + + +def self_cal(visibilities, expected_sources, baselines, frequency, iterations=10): + model = estimate_model_visibilities(expected_sources, visibilities, baselines, frequency) + for i in range(iterations): + model, residual = compute_calibrated_model(visibilities, model, maxiter=50) + logging.debug('residual after iteration %d: %s', i, residual) + return model + + def simulate_sky_source(lmn_coord: np.array, baselines: np.array, freq: float): """ Simulate visibilities for a sky source diff --git a/lofarimaging/singlestationutil.py b/lofarimaging/singlestationutil.py index 0b57b43..995f4b2 100644 --- a/lofarimaging/singlestationutil.py +++ b/lofarimaging/singlestationutil.py @@ -1,36 +1,32 @@ """Functions for working with LOFAR single station data""" -import os import datetime +import os from typing import List, Dict, Tuple, Union -import numpy as np -from packaging import version -import tqdm +import astropy.units as u import h5py - -import matplotlib.pyplot as plt +import lofarantpos +import lofargeotiff import matplotlib.animation -from matplotlib.ticker import FormatStrFormatter -from matplotlib import cm -from matplotlib.figure import Figure -from matplotlib.colors import ListedColormap, Normalize -from matplotlib.patches import Circle import matplotlib.axes as maxes -from mpl_toolkits.axes_grid1 import make_axes_locatable - +import matplotlib.pyplot as plt +import numpy as np +import tqdm from astropy.coordinates import SkyCoord, GCRS, EarthLocation, AltAz, get_sun, get_moon -import astropy.units as u from astropy.time import Time - -import lofargeotiff from lofarantpos.db import LofarAntennaDatabase -import lofarantpos +from matplotlib import cm +from matplotlib.colors import ListedColormap, Normalize +from matplotlib.figure import Figure +from matplotlib.patches import Circle +from matplotlib.ticker import FormatStrFormatter +from mpl_toolkits.axes_grid1 import make_axes_locatable +from packaging import version -from .maputil import get_map, make_leaflet_map -from .lofarimaging import nearfield_imager, sky_imager, skycoord_to_lmn, subtract_sources from .hdf5util import write_hdf5 - +from .lofarimaging import nearfield_imager, sky_imager, skycoord_to_lmn, subtract_sources +from .maputil import get_map, make_leaflet_map __all__ = ["sb_from_freq", "freq_from_sb", "find_caltable", "read_caltable", "rcus_in_station", "read_acm_cube", "get_station_pqr", "get_station_type", @@ -331,7 +327,7 @@ def get_station_pqr(station_name: str, rcu_mode: Union[str, int], db): 'intl': GENERIC_INT_201512, 'remote': GENERIC_REMOTE_201512, 'core': GENERIC_CORE_201512 } selected_dipoles = selected_dipole_config[station_type] + \ - np.arange(len(selected_dipole_config[station_type])) * 16 + np.arange(len(selected_dipole_config[station_type])) * 16 station_pqr = db.hba_dipole_pqr(full_station_name)[selected_dipoles] else: raise RuntimeError("Station name did not contain LBA or HBA, could not load antenna positions") @@ -340,7 +336,7 @@ def get_station_pqr(station_name: str, rcu_mode: Union[str, int], db): def make_ground_plot(image: np.ndarray, background_map: np.ndarray, extent: List[float], title: str = "Ground plot", - subtitle: str = "", opacity: float = 0.6, fig: Figure = None, draw_contours: bool = True, **kwargs) \ + subtitle: str = "", opacity: float = 0.6, fig: Figure = None, draw_contours: bool = True, **kwargs) \ -> Tuple[Figure, np.ndarray]: """ Make a ground plot of an array with data @@ -585,7 +581,6 @@ def make_xst_plots(xst_data: np.ndarray, opacity: Opacity for map overlay. Defaults to 0.6. hdf5_filename: Filename where hdf5 results can be written. Defaults to outputpath + '/results.h5' outputpath: Directory where results can be saved. Defaults to 'results' - subtract: List of sources to subtract. Defaults to None Returns: @@ -662,14 +657,14 @@ def make_xst_plots(xst_data: np.ndarray, marked_bodies = { 'Cas A': SkyCoord(ra=350.85 * u.deg, dec=58.815 * u.deg), 'Cyg A': SkyCoord(ra=299.86815191 * u.deg, dec=40.73391574 * u.deg), - 'Per A': SkyCoord(ra=49.95066567*u.deg, dec=41.51169838 * u.deg), - 'Her A': SkyCoord(ra=252.78343333*u.deg, dec=4.99303056*u.deg), - 'Cen A': SkyCoord(ra=201.36506288*u.deg, dec=-43.01911267*u.deg), - 'Vir A': SkyCoord(ra=187.70593076*u.deg, dec=12.39112329*u.deg), - '3C295': SkyCoord(ra=212.83527917*u.deg, dec=52.20264444*u.deg), + 'Per A': SkyCoord(ra=49.95066567 * u.deg, dec=41.51169838 * u.deg), + 'Her A': SkyCoord(ra=252.78343333 * u.deg, dec=4.99303056 * u.deg), + 'Cen A': SkyCoord(ra=201.36506288 * u.deg, dec=-43.01911267 * u.deg), + 'Vir A': SkyCoord(ra=187.70593076 * u.deg, dec=12.39112329 * u.deg), + '3C295': SkyCoord(ra=212.83527917 * u.deg, dec=52.20264444 * u.deg), 'Moon': get_moon(obstime_astropy, location=station_earthlocation).transform_to(GCRS), 'Sun': get_sun(obstime_astropy), - '3C196': SkyCoord(ra=123.40023371*u.deg, dec=48.21739888*u.deg) + '3C196': SkyCoord(ra=123.40023371 * u.deg, dec=48.21739888 * u.deg) } marked_bodies_lmn = {} @@ -752,7 +747,7 @@ def make_xst_plots(xst_data: np.ndarray, "pixels_per_metre": pixels_per_metre} tags.update(calibration_info) lon_min, lon_max, lat_min, lat_max = extent_lonlat - lofargeotiff.write_geotiff(ground_img[::-1,:], os.path.join(outputpath, f"{fname}_nearfield_calibrated.tiff"), + lofargeotiff.write_geotiff(ground_img[::-1, :], os.path.join(outputpath, f"{fname}_nearfield_calibrated.tiff"), (lon_min, lat_max), (lon_max, lat_min), as_pqr=False, stationname=station_name, obsdate=obstime, tags=tags) @@ -769,7 +764,7 @@ def make_sky_movie(moviefilename: str, h5file: h5py.File, obsnums: List[str], vm """ Make movie of a list of observations """ - fig = plt.figure(figsize=(10,10)) + fig = plt.figure(figsize=(10, 10)) for obsnum in tqdm.tqdm(obsnums): obs_h5 = h5file[obsnum] skydata_h5 = obs_h5["sky_img"] @@ -787,7 +782,28 @@ def make_sky_movie(moviefilename: str, h5file: h5py.File, obsnums: List[str], vm # Thanks to Maaijke Mevius for making this animation work! ims = fig.get_children()[1:] - ims = [ims[i:i+2] for i in range(0, len(ims), 2)] + ims = [ims[i:i + 2] for i in range(0, len(ims), 2)] ani = matplotlib.animation.ArtistAnimation(fig, ims, interval=30, blit=False, repeat_delay=1000) writer = matplotlib.animation.writers['ffmpeg'](fps=5, bitrate=800) ani.save(moviefilename, writer=writer, dpi=fig.dpi) + + +def compute_baselines(station_name, rcu_mode): + # Setup the database + db = LofarAntennaDatabase() + + station_pqr = get_station_pqr(station_name, rcu_mode, db) + + station_name = get_full_station_name(station_name, rcu_mode) + + # Rotate station_pqr to a north-oriented xyz frame, where y points North, in a plane through the station. + rotation = db.rotation_from_north(station_name) + + pqr_to_xyz = np.array([[np.cos(-rotation), -np.sin(-rotation), 0], + [np.sin(-rotation), np.cos(-rotation), 0], + [0, 0, 1]]) + + station_xyz = (pqr_to_xyz @ station_pqr.T).T + + baselines = station_xyz[:, np.newaxis, :] - station_xyz[np.newaxis, :, :] + return baselines diff --git a/self_cal_demo.py b/self_cal_demo.py new file mode 100644 index 0000000..24082ac --- /dev/null +++ b/self_cal_demo.py @@ -0,0 +1,66 @@ +import h5py +import numpy +import lofarimaging.hdf5util as h5utils +from lofarimaging.singlestationutil import compute_baselines, make_sky_plot +from lofarimaging.lofarimaging import compute_pointing_matrix, compute_calibrated_model, sky_imager,\ + self_cal, estimate_model_visibilities +import matplotlib.pyplot as plt + +test_data_path = '/home/mmancini/Documents/Projects/SingleStationImager/test_dataset.h5' +test_station = 'CS103' +test_data = h5py.File(test_data_path, 'r') +obs_names = h5utils.get_obsnums(test_data, station_name=test_station) +an_obs = test_data[obs_names[10]] + +a_frequency = an_obs.attrs['frequency'] +a_source_lmn = an_obs.attrs['source_lmn'] +a_source_name = an_obs.attrs['source_names'] +a_rcu_mode = an_obs.attrs['rcu_mode'] +a_dataset = an_obs['calibrated_data'] + +a_dataset_xx = a_dataset[::2, ::2] +a_dataset_yy = a_dataset[1::2, 1::2] + +a_dataset_i = a_dataset_xx + a_dataset_yy +a_dataset_q = a_dataset_xx - a_dataset_yy +a_dataset_u = 2 * (a_dataset_xx * a_dataset_yy.conj()).real +a_dataset_v = -2 * (a_dataset_xx * a_dataset_yy.conj()).imag + +baselines = compute_baselines(test_station, a_rcu_mode) +sky_image = sky_imager(a_dataset_i, baselines, a_frequency, 100, 100) + +model = estimate_model_visibilities(a_source_lmn, a_dataset_i, baselines, a_frequency) + +model_image = sky_imager(model, baselines, a_frequency, 100, 100) + +pointing_matrix = compute_pointing_matrix(a_source_lmn, baselines, a_frequency) +calibrated_model, _ = compute_calibrated_model(a_dataset_i, model, maxiter=50) +self_cal_model = self_cal(a_dataset_i, a_source_lmn, baselines, a_frequency) + +sky_without_calibrated_model = sky_imager(a_dataset_i - self_cal_model, baselines, a_frequency, 100, 100) +calibrated_model_image = sky_imager(calibrated_model, baselines, a_frequency, 100, 100) +sky_without_model = sky_imager(a_dataset_i - model, baselines, a_frequency, 100, 100) + + +f, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4) +min, max = numpy.min(sky_image), numpy.max(sky_image) + +ax1.scatter(a_source_lmn[:, 0], a_source_lmn[:, 1], marker='o', color='r', s=100, facecolors='none') +im = ax1.imshow(sky_image[::-1,:], vmin=min, vmax=max, extent=(1, -1, -1, 1)) +ax1.set_title('initial sky') + +ax2.scatter(a_source_lmn[:, 0], a_source_lmn[:, 1], marker='o', color='r', s=100, facecolors='none') +ax2.imshow(sky_without_model[::-1, :], extent=(1, -1, -1, 1), vmin=min, vmax=max) +ax2.set_title('sky without model') + +ax3.scatter(a_source_lmn[:, 0], a_source_lmn[:, 1], marker='o', color='r', s=100, facecolors='none') +ax3.imshow(calibrated_model_image[::-1, :], extent=(1, -1, -1, 1)) +ax3.set_title('calibrated model image') + +ax4.imshow(sky_without_calibrated_model[::-1, :], extent=(1, -1, -1, 1), vmin=min, vmax=max) +to_plot = {name: pos for pos, name in zip(a_source_lmn, a_source_name)} +ax4.set_title('sky without calibrated model') + +make_sky_plot(sky_image, to_plot) + +plt.show()