Skip to content
Snippets Groups Projects
Commit ebd5874e authored by Jakob Maljaars's avatar Jakob Maljaars
Browse files

Document

parent bc7abeac
No related branches found
No related tags found
1 merge request!58Refactor idg cal
Pipeline #8490 passed
Stage: prepare
Stage: build
Stage: integration_and_deploy
Show whitespace changes
Inline Side-by-side
idg-cal/basisfunctions.py
+ 69
− 40
View file @ ebd5874e
# 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()
idg-cal/idgcaldpstep.py
+ 80
− 36
View file @ ebd5874e
......@@ -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
idg-cal/unit_tests/test_basisfunctions.py
+ 6
− 4
View file @ ebd5874e
#!/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))),
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment