From ebd5874e68e84339a3f7a47c22ca076eefba4a81 Mon Sep 17 00:00:00 2001 From: Jakob Maljaars <jakob.maljaars@stcorp.nl> Date: Tue, 26 Jan 2021 15:54:51 +0100 Subject: [PATCH] Document --- idg-cal/basisfunctions.py | 109 ++++++++++++-------- idg-cal/idgcaldpstep.py | 116 +++++++++++++++------- idg-cal/unit_tests/test_basisfunctions.py | 10 +- 3 files changed, 155 insertions(+), 80 deletions(-) diff --git a/idg-cal/basisfunctions.py b/idg-cal/basisfunctions.py index a44dd61fc..8ff356666 100644 --- a/idg-cal/basisfunctions.py +++ b/idg-cal/basisfunctions.py @@ -1,9 +1,34 @@ +# Copyright (C) 2020 ASTRON (Netherlands Institute for Radio Astronomy) +# SPDX-License-Identifier: GPL-3.0-or-later + import numpy as np from numpy.polynomial.polynomial import polyval as np_polyval class LagrangePolynomial: + """ + Class implementing functionality for (Lagrange)Polynomials + """ + def __init__(self, order=None, nr_coeffs=None): + """ + Initialize LagrangePolynomial, by either specifying + + Parameters + ---------- + order : int, optional + Polynomial order. If not specified, it will be + derived from the number of coefficients + nr_coeffs : int, optional + Number of coefficients. If not specified, it will be + derived from the polynomial order. + + Raises + ------ + ValueError + If neither the order nor the nr_coeffs are specified, + or if both are specified. + """ if (order is None and nr_coeffs is None) or ( order is not None and nr_coeffs is not None ): @@ -22,17 +47,17 @@ class LagrangePolynomial: Parameters ---------- - x : [type] - [description] - y : [type] - [description] - coeffs : [type] - [description] + x : float, np.ndarray + x-coordinate(s), can be scalar or 1D vector + y : float, np.ndarray + y-coordinate(s), can be scalar or 1D vector + coeffs : np.ndarray + Expansion coefficients Returns ------- - [type] - [description] + np.ndarray + Array storing the polynomial evaluations at x,y """ # Evaluate polynomial using Horner's method assert coeffs.size == self.nr_coeffs @@ -53,24 +78,23 @@ class LagrangePolynomial: def expand_basis(self, x, y): """ - Expand the basis on the give coordinates + Expand the basis for the given coordinates. Accepted input is: + - x and y both scalar + - x and y both vector (which will be expanded on a grid) Parameters ---------- - x : scalar, np.ndarray - [description] - y : scalar, np.ndarray - [description] + x : float, np.ndarray + x coordinate(s) + y : float, np.ndarray + y coordinate(s) Returns ------- - [type] - [description] - - Raises - ------ - TypeError - [description] + np.ndarray + numpy array with expanded basis functions. Return size will be: + - if x and y both scalar (nr_coeffs, 1) array + - if x and y both vector (nr_coeffs, x.len, y.len) """ if not isinstance(x, np.ndarray) and not isinstance(y, np.ndarray): @@ -118,6 +142,31 @@ class LagrangePolynomial: @staticmethod def get_indices_right_diagonal(order, diag_nr): + """ + Get indices in corresponding to right diagonal number in Pascal's triangle. + Take the triangle for a second order polynomial + + :: + 0 + 1 2 + 3 4 5 + + Diagonal 2 will return [3], diagonal 1 will return [1, 4] and + diagonal 0 will return [5, 2, 0] + + Parameters + ---------- + order : int + Polynomial order + diag_nr : int + Index of the diagonal. Should be <= order + + Returns + ------- + np.ndarray + Numpy array with indices + """ + if diag_nr > order: raise ValueError("Diagonal number should be smaller than polynomial order") @@ -129,23 +178,3 @@ class LagrangePolynomial: for i in range(diag_nr + 1, order + 1): indices.append(indices[-1] + i + 1) return np.array(indices) - - -# TODO: make unit tests -def main(): - poly = LagrangePolynomial(order=1) - # result = poly.evaluate(1, 2, coeffs=np.array([1, 2, 3])) - x = 1 - y = np.array([1, 2]) - - # X, Y = np.meshgrid(x, y) - # print(X) - # print(Y) - result = poly.evaluate(x, y, coeffs=np.array([1, 2, 3])) - # result = poly.evaluate(2, 3, coeffs=np.array([1, 2, 3, 4, 5, 6])) - print(result) - # print(np_polyval()) - - -if __name__ == "__main__": - main() diff --git a/idg-cal/idgcaldpstep.py b/idg-cal/idgcaldpstep.py index a2b9fe80d..be00669f0 100644 --- a/idg-cal/idgcaldpstep.py +++ b/idg-cal/idgcaldpstep.py @@ -298,34 +298,15 @@ class IDGCalDPStep(dppp.DPStep): self.initialize() # Concatenate accumulated data and display just the shapes - visibilities = np.stack( - [np.array(dpbuffer.get_data(), copy=False) for dpbuffer in self.dpbuffers], - axis=1, - ) - visibilities = visibilities[self.auto_corr_mask, :, :] - flags = np.stack( - [np.array(dpbuffer.get_flags(), copy=False) for dpbuffer in self.dpbuffers], - axis=1, - ) - flags = flags[self.auto_corr_mask, :, :] - weights = np.stack( - [ - np.array(dpbuffer.get_weights(), copy=False) - for dpbuffer in self.dpbuffers - ], - axis=1, - ) - weights = weights[self.auto_corr_mask, :, :] - uvw_ = np.stack( - [np.array(dpbuffer.get_uvw(), copy=False) for dpbuffer in self.dpbuffers], - axis=1, - ) + visibilities = self.__extract_buffer("visibilities") + flags = self.__extract_buffer("flags") + weights = self.__extract_buffer("weights") + uvw_ = self.__extract_buffer("uvw", apply_autocorr_mask=False) uvw = np.zeros(shape=(self.nr_baselines, self.ampl_interval), dtype=idg.uvwtype) print(f" shape uvw_: {uvw_.shape}") print(self.auto_corr_mask.shape) - print(uvw_[self.auto_corr_mask, :, 0].shape) uvw["u"] = uvw_[self.auto_corr_mask, :, 0] uvw["v"] = -uvw_[self.auto_corr_mask, :, 1] uvw["w"] = -uvw_[self.auto_corr_mask, :, 2] @@ -376,6 +357,9 @@ class IDGCalDPStep(dppp.DPStep): aterm_offsets = get_aterm_offsets(self.nr_phase_updates, self.ampl_interval) + # parameters has shape (nr_stations, nr_coeffs) for amplitude and (nr_stations, nr_phase_updates, nr_coeffs) + # for amplitude + # amplitude/phase basis (Bampl/Bphase) (nr_coeffs, subgridsize, subgridzise, nr_polarizations) aterm_ampl = np.tensordot( parameters[:, : self.ampl_poly.nr_coeffs], self.Bampl, axes=((1,), (0,)) ) @@ -390,6 +374,7 @@ class IDGCalDPStep(dppp.DPStep): ) ) + # aterms will have shape (nr_phase_updates, nr_stations, subgrid size, subgrid size, nr polarizations) aterms = np.ascontiguousarray( (aterm_phase.transpose((1, 0, 2, 3, 4)) * aterm_ampl).astype( idg.idgtypes.atermtype @@ -398,16 +383,13 @@ class IDGCalDPStep(dppp.DPStep): nr_iterations = 0 converged = False + previous_residual = 0.0 max_dx = 0.0 timer = -time.time() - timer0 = 0 timer1 = 0 - - previous_residual = 0.0 - while True: nr_iterations += 1 print(f"iteration nr {nr_iterations} ") @@ -416,7 +398,7 @@ class IDGCalDPStep(dppp.DPStep): norm_dx = 0.0 residual_sum = 0.0 for i in range(self.nr_stations): - print(f" {i}") + print(f" Station {i}") timer1 -= time.time() # Predict visibilities for current solution @@ -449,6 +431,7 @@ class IDGCalDPStep(dppp.DPStep): ) ) + # new-axis is introduced at "stations" axis aterm_derivatives_ampl = ( aterm_phase[:, np.newaxis, :, :, :] * self.Bampl[np.newaxis, :, :, :, :] @@ -474,9 +457,7 @@ class IDGCalDPStep(dppp.DPStep): ) timer0 += time.time() - residual0 = residual[0] - residual_sum += residual[0] gradient = np.concatenate( @@ -530,10 +511,8 @@ class IDGCalDPStep(dppp.DPStep): i_max = i p0 = parameters[i].copy() - parameters[i] += self.solver_update_gain * dx - # TODO: probably no need to repeat when writing to H5Parm aterm_ampl = np.repeat( np.tensordot( parameters[i, : self.ampl_poly.nr_coeffs], @@ -594,7 +573,7 @@ class IDGCalDPStep(dppp.DPStep): amplitude_coefficients = parameters_polynomial[ :, : self.ampl_poly.nr_coeffs ].reshape(self.nr_stations, 1, self.ampl_poly.nr_coeffs) - # phase parameters: reshaped into (nr_stations, nr_timeslots, nr_parameters_phase) array + # phase parameters: reshaped into (nr_stations, nr_phase_updates, nr_parameters_phase) array phase_coefficients = parameters_polynomial[ :, self.ampl_poly.nr_coeffs : : ].reshape(self.nr_stations, self.nr_phase_updates, self.phase_poly.nr_coeffs) @@ -607,9 +586,47 @@ class IDGCalDPStep(dppp.DPStep): self.h5writer.fill_solution_table( "phase_coefficients", phase_coefficients, offset_phase ) - self.count_process_buffer_calls += 1 + def __extract_buffer(self, name, apply_autocorr_mask=True): + """ + Extract buffer from buffered data. + + Parameters + ---------- + name : str + Should be any of ("visibilities", "weights", "flags", "uvw") + apply_autocorr_mask : bool, optional + Remove autocorrelation from returned result? Defaults to True + + Returns + ------- + np.ndarray + """ + + if name == "visibilities": + result = [ + np.array(dpbuffer.get_data(), copy=False) for dpbuffer in self.dpbuffers + ] + elif name == "flags": + result = [ + np.array(dpbuffer.get_flags(), copy=False) + for dpbuffer in self.dpbuffers + ] + elif name == "weights": + result = [ + np.array(dpbuffer.get_weights(), copy=False) + for dpbuffer in self.dpbuffers + ] + elif name == "uvw": + result = [ + np.array(dpbuffer.get_uvw(), copy=False) for dpbuffer in self.dpbuffers + ] + else: + raise ValueError("Name not recognized") + result = np.stack(result, axis=1) + return result[self.auto_corr_mask, :, :] if apply_autocorr_mask else result + def init_h5parm_solution_table( h5parm_object, @@ -622,6 +639,35 @@ def init_h5parm_solution_table( basisfunction_type="lagrange", history="", ): + """ + Initialize h5parm solution table + + Parameters + ---------- + h5parm_object : idg.h5parmwriter.H5ParmWriter + h5parm object + soltab_type : str + Any of ("amplitude", "phase") + axes_info : dict + Dict containing axes info (name, length) + antenna_names : np.ndarray + Array of strings containing antenna names + time_array : np.ndarray + Array of times + image_size : float + Pixel size + subgrid_size : int + Subgrid size, used in IDG + basisfunction_type : str, optional + Which basis function was used? Defaults to "lagrange" + history : str, optional + History attribute, by default "" + + Returns + ------- + idg.h5parmwriter.H5ParmWriter + Extended H5ParmWriter object + """ soltab_info = {"amplitude": "amplitude_coefficients", "phase": "phase_coefficients"} assert soltab_type in soltab_info.keys() @@ -834,8 +880,6 @@ def idgwindow(N, W, padding, offset=0.5, l_range=None): def get_aterm_offsets(nr_timeslots, nr_time): aterm_offsets = np.zeros((nr_timeslots + 1), dtype=idg.idgtypes.atermoffsettype) - for i in range(nr_timeslots + 1): aterm_offsets[i] = i * (nr_time // nr_timeslots) - return aterm_offsets diff --git a/idg-cal/unit_tests/test_basisfunctions.py b/idg-cal/unit_tests/test_basisfunctions.py index 065f0d832..de022cbb7 100644 --- a/idg-cal/unit_tests/test_basisfunctions.py +++ b/idg-cal/unit_tests/test_basisfunctions.py @@ -1,3 +1,7 @@ +#!/usr/bin/env python3 +# Copyright (C) 2020 ASTRON (Netherlands Institute for Radio Astronomy) +# SPDX-License-Identifier: GPL-3.0-or-later3 + import sys sys.path.append("..") @@ -5,8 +9,6 @@ from basisfunctions import LagrangePolynomial import numpy as np import pytest -# import os - def explicit_evaluation(x, y, order, coeffs): sol = 0.0 @@ -26,9 +28,9 @@ def explicit_evaluation(x, y, order, coeffs): (1.0, 2.0, np.array(range(1, 4))), (1.0, 2.0, np.array(range(1, 7))), (1.0, 2.0, np.array(range(1, 11))), - # # x vector, y scalar (only third order polynomial) + # x vector, y scalar (only third order polynomial) (np.array(range(-2, 10)), 2.0, np.array(range(1, 11))), - # # x scalar, y vector (only third order polynomial) + # x scalar, y vector (only third order polynomial) (2.0, np.array(range(-2, 10)), np.array(range(1, 11))), # both vector (outcome is 2d matrix) (np.array(range(-2, 10)), np.array(range(-2, 10)), np.array(range(1, 11))), -- GitLab