diff --git a/idg-cal/basisfunctions.py b/idg-cal/basisfunctions.py
index a44dd61fc0268c68310478b22891d7f0b254b960..8ff356666095b4369ed870c850a117642feeafd9 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 a2b9fe80d61b3aaa8141c74e1ed6bcd18d0ab381..be00669f0e7b27108e0277ddf073b4e1a2093dde 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 065f0d8324164fdacbd0d0becbd351df8e948839..de022cbb79684f20284bbaeed7d22b44830848f0 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))),