diff --git a/.gitattributes b/.gitattributes
index 540936775fc05646fb55638ba57dbde45fc36b4a..65929653f8ee137bd0b4287cdb2057c1db64e44f 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -40,13 +40,16 @@ CAL/CalibrationProcessing/lib/__init__.py -text
CAL/CalibrationProcessing/lib/process_holography_dataset.py -text
CAL/CalibrationProcessing/lib/processing/__init__.py -text
CAL/CalibrationProcessing/lib/processing/averaging.py -text
-CAL/CalibrationProcessing/lib/processing/matrices.py -text
+CAL/CalibrationProcessing/lib/processing/normalize.py -text
+CAL/CalibrationProcessing/lib/processing/solver.py -text
CAL/CalibrationProcessing/test/CMakeLists.txt -text
+CAL/CalibrationProcessing/test/__init__.py -text
+CAL/CalibrationProcessing/test/processing/CMakeLists.txt -text
+CAL/CalibrationProcessing/test/processing/t_processing.py -text
+CAL/CalibrationProcessing/test/processing/t_processing.run -text
+CAL/CalibrationProcessing/test/processing/t_processing.sh -text
CAL/CalibrationProcessing/test/profiling_processing.py -text
-CAL/CalibrationProcessing/test/t_processing.py -text
-CAL/CalibrationProcessing/test/t_processing.run -text
-CAL/CalibrationProcessing/test/t_processing.sh -text
-CAL/CalibrationProcessing/test/test_utils.py -text
+CAL/CalibrationProcessing/test/utils.py -text
CEP/Calibration/BBSControl/include/BBSControl/CommandHandlerEstimator.h -text
CEP/Calibration/BBSControl/include/BBSControl/CommandHandlerReducer.h -text
CEP/Calibration/BBSControl/include/BBSControl/OptionParser.h -text
diff --git a/CAL/CalibrationCommon/lib/datacontainers/__init__.py b/CAL/CalibrationCommon/lib/datacontainers/__init__.py
index 33d69fa5d3fde3cbfe6a3f2510ebe61d266ae825..7418ea83e6c5a1715534210508bd2b282fdf9f90 100644
--- a/CAL/CalibrationCommon/lib/datacontainers/__init__.py
+++ b/CAL/CalibrationCommon/lib/datacontainers/__init__.py
@@ -1,8 +1,8 @@
+__all__ = ['HolographyDataset', 'HolographyObservation', 'HolographySpecification', 'HolographyMeasurementSet']
+
from .holography_specification import HolographySpecification
from .holography_observation import HolographyObservation
from .holography_measurementset import HolographyMeasurementSet
from .holography_dataset import HolographyDataset
-
-__all__ = ['HolographyDataset', 'HolographyObservation', 'HolographySpecification', 'HolographyMeasurementSet']
\ No newline at end of file
diff --git a/CAL/CalibrationProcessing/lib/CMakeLists.txt b/CAL/CalibrationProcessing/lib/CMakeLists.txt
index 38c673e25925f17c0af28a11d64af92dba1972c3..4434b356d113279e9ed80303475315ad096851e8 100644
--- a/CAL/CalibrationProcessing/lib/CMakeLists.txt
+++ b/CAL/CalibrationProcessing/lib/CMakeLists.txt
@@ -20,8 +20,9 @@
python_install(
__init__.py
processing/averaging.py
- processing/matrices.py
+ processing/normalize.py
+ processing/solver.py
processing/__init__.py
DESTINATION lofar/calibration)
-lofar_add_bin_script(process_holography_dataset.py RENAME process_holography_dataset)
\ No newline at end of file
+lofar_add_bin_script(process_holography_dataset.py RENAME process_holography_dataset)
diff --git a/CAL/CalibrationProcessing/lib/processing/__init__.py b/CAL/CalibrationProcessing/lib/processing/__init__.py
index a26b7d7eb0ddb5d3ebae74aacc521aea1e25953b..90a414f8660ea4a6a269820a8fb54ee89bff8572 100644
--- a/CAL/CalibrationProcessing/lib/processing/__init__.py
+++ b/CAL/CalibrationProcessing/lib/processing/__init__.py
@@ -1,2 +1,2 @@
from .averaging import *
-from .matrices import *
\ No newline at end of file
+from .normalize import *
\ No newline at end of file
diff --git a/CAL/CalibrationProcessing/lib/processing/averaging.py b/CAL/CalibrationProcessing/lib/processing/averaging.py
index 11d031c92c4dd6a7010710aca76867df6f21f198..56f518a660750b827996e6f12613ac68cf9f7127 100644
--- a/CAL/CalibrationProcessing/lib/processing/averaging.py
+++ b/CAL/CalibrationProcessing/lib/processing/averaging.py
@@ -1,8 +1,9 @@
import logging
import numpy
-
from lofar.calibration.common.datacontainers import HolographyDataset
+from lofar.calibration.common.datacontainers.holography_dataset_definitions \
+ import HDS_data_sample_type
logger = logging.getLogger(__name__)
@@ -178,7 +179,7 @@ def average_datatable_by_time_interval(data_table_in, time_interval):
def _average_dataset_by_function(dataset, function, *parameters):
out_data = {}
- for reference_station, data_per_reference_station in dataset.data.items():
+ for reference_station, data_per_reference_station in dataset.items():
if reference_station not in out_data:
out_data[reference_station] = dict()
@@ -234,16 +235,17 @@ def average_dataset_by_sample(input_data_set, window_size):
logger.error(text)
raise ValueError(text)
else:
- text = "Cannot average data by samples! The input data is of the invalid type %s." % (type(input_data_set))
- logger.error(text)
- raise ValueError(text)
+ text = "Cannot average data by samples! The input data is of the invalid type %s." % (
+ type(input_data_set))
+ logger.error(text)
+ raise ValueError(text)
if isinstance(int, window_size) is False or window_size < 1:
text = "Cannot average data by samples! The window size must be an integer value bigger than 0."
logger.error(text)
raise ValueError(text)
- out_data = _average_dataset_by_function(dataset, average_datatable_by_sample, window_size)
+ out_data = _average_dataset_by_function(dataset.data, average_datatable_by_sample, window_size)
if dataset.derived_data is None:
dataset.derived_data = dict()
@@ -279,16 +281,17 @@ def average_dataset_by_time(input_data_set, time_interval):
logger.error(text)
raise ValueError(text)
else:
- text = "Cannot average data by time! The input data is of the invalid type %s." % (type(input_data_set))
- logger.error(text)
- raise ValueError(text)
+ text = "Cannot average data by time! The input data is of the invalid type %s." % (
+ type(input_data_set))
+ logger.error(text)
+ raise ValueError(text)
if isinstance(float, time_interval) is False or time_interval < 0:
text = "Cannot average data by samples! The time interval must be positive."
logger.error(text)
raise ValueError(text)
- out_data = _average_dataset_by_function(dataset, average_datatable_by_time_interval,
+ out_data = _average_dataset_by_function(dataset.data, average_datatable_by_time_interval,
time_interval)
if dataset.derived_data is None:
dataset.derived_data = dict()
@@ -298,6 +301,12 @@ def average_dataset_by_time(input_data_set, time_interval):
return dataset
+def average_data_by_sample(input_data_table, samples):
+ out_data = _average_dataset_by_function(input_data_table, average_datatable_by_sample,
+ samples)
+ return out_data
+
+
def _compute_size_of_new_array(array_size, window_size):
"""
Round up the array size given a window size.
@@ -315,3 +324,87 @@ def _compute_size_of_new_array(array_size, window_size):
# add one if the len of the data_array is not divisible by the sample_width
new_array_size += 1 if array_size % window_size else 0
return new_array_size
+
+
+def _weigthed_average_data_per_polarization(input_dataset, output_dataset,
+ real_weights, imag_weights, polarization):
+ mask = input_dataset['flag'] == False
+ try:
+ output_dataset[polarization] = numpy.average(
+ input_dataset[polarization][mask].real,
+ weights=real_weights[polarization][mask]) \
+ + 1.j * numpy.average(
+ input_dataset[polarization][mask].imag,
+ weights=imag_weights[polarization][mask])
+ except ZeroDivisionError as e:
+ logger.warning('error in weighted average : %s', e)
+ output_dataset['flag'] = True
+
+def weighted_average_dataset_per_station(dataset, input_data_table):
+ """
+
+ :param dataset:
+ :type dataset: HolographyDataset
+ :param input_data_table:
+ :return:
+ """
+ stations = dataset.reference_stations
+ n_stations = len(stations)
+ frequencies = dataset.frequencies
+ beams = dataset.beamlets
+
+ result_ALL = dict()
+ result = dict(all=result_ALL)
+ for frequency in frequencies:
+ result_per_frequency = dict()
+ result_ALL[frequency] = result_per_frequency
+
+ for beam in beams:
+ beam_str = str(beam)
+ result_per_beam = numpy.zeros(1, dtype=HDS_data_sample_type)
+
+
+ average_per_beam_station = numpy.zeros(n_stations, dtype=HDS_data_sample_type)
+ std_real_per_beam_station = numpy.zeros(n_stations, dtype=HDS_data_sample_type)
+ std_imag_per_beam_station = numpy.zeros(n_stations, dtype=HDS_data_sample_type)
+
+ for i, station in enumerate(stations):
+ data_per_station_frequency_beam = \
+ input_data_table[station][str(frequency)][beam_str]
+ average_per_beam_station['XX'][i] = data_per_station_frequency_beam['mean']['XX']
+ average_per_beam_station['XY'][i] = data_per_station_frequency_beam['mean']['XY']
+ average_per_beam_station['YX'][i] = data_per_station_frequency_beam['mean']['YX']
+ average_per_beam_station['YY'][i] = data_per_station_frequency_beam['mean']['YY']
+
+ average_per_beam_station['l'][i] = data_per_station_frequency_beam['mean']['l']
+ average_per_beam_station['m'][i] = data_per_station_frequency_beam['mean']['m']
+
+ std_real_per_beam_station['XX'][i] = \
+ data_per_station_frequency_beam['std']['XX_real']
+ std_real_per_beam_station['XY'][i] = \
+ data_per_station_frequency_beam['std']['XY_real']
+ std_real_per_beam_station['YX'][i] = \
+ data_per_station_frequency_beam['std']['YX_real']
+ std_real_per_beam_station['YY'][i] = \
+ data_per_station_frequency_beam['std']['YY_real']
+
+ std_imag_per_beam_station['XX'][i] = \
+ data_per_station_frequency_beam['std']['XX_imag']
+ std_imag_per_beam_station['XY'][i] = \
+ data_per_station_frequency_beam['std']['XY_imag']
+ std_imag_per_beam_station['YX'][i] = \
+ data_per_station_frequency_beam['std']['YX_imag']
+ std_imag_per_beam_station['YY'][i] = \
+ data_per_station_frequency_beam['std']['YY_imag']
+
+ for polarization in ['XX', 'XY', 'YX', 'YY']:
+ _weigthed_average_data_per_polarization(average_per_beam_station,
+ result_per_beam,
+ std_real_per_beam_station,
+ std_imag_per_beam_station,
+ polarization)
+ result_per_beam['l'] = numpy.average(average_per_beam_station['l'])
+ result_per_beam['m'] = numpy.average(average_per_beam_station['m'])
+
+ result_per_frequency[beam_str] = dict(mean=result_per_beam)
+ return result
diff --git a/CAL/CalibrationProcessing/lib/processing/matrices.py b/CAL/CalibrationProcessing/lib/processing/matrices.py
deleted file mode 100644
index 232678327b5c67240c464230cb97d96c8c26104d..0000000000000000000000000000000000000000
--- a/CAL/CalibrationProcessing/lib/processing/matrices.py
+++ /dev/null
@@ -1,29 +0,0 @@
-import logging
-
-import numpy
-
-logger = logging.getLogger(__name__)
-
-
-def extract_crosscorrelation_matrices_from_HDS_datatable(datatable):
- new_shape = (2, 2, datatable.shape[0])
- new_array = numpy.zeros(shape=new_shape, dtype=numpy.complex)
- new_array[0, 0, :] = datatable['XX']
- new_array[0, 1, :] = datatable['XY']
- new_array[1, 0, :] = datatable['YX']
- new_array[1, 1, :] = datatable['YY']
-
- return new_array
-
-
-def invert_crosscorrelation_matrices(cross_correlation_matrices):
- """
- Invert the matrices along the last axis
- :param cross_correlation_matrices: matrices to invert
- :type cross_correlation_matrices: numpy.ndarray
- :return:
- """
- assert cross_correlation_matrices.ndim >= 3
-
- return numpy.linalg.inv(cross_correlation_matrices)
-
diff --git a/CAL/CalibrationProcessing/lib/processing/normalize.py b/CAL/CalibrationProcessing/lib/processing/normalize.py
new file mode 100644
index 0000000000000000000000000000000000000000..ba43ad48b83ae12937e14e114a2869e36d32209d
--- /dev/null
+++ b/CAL/CalibrationProcessing/lib/processing/normalize.py
@@ -0,0 +1,135 @@
+import logging
+
+import numpy
+
+logger = logging.getLogger(__name__)
+
+__POLARIZATION_TO_INDEX = {
+ 'XX': (0,0),
+ 'XY': (0,1),
+ 'YX': (1,0),
+ 'YY': (1,1)
+}
+
+
+def normalize_beams_by_central_one(input_data, central_beamlet_number):
+ """
+
+ :param central_beamlet_number: number of the central beamleat
+ :type central_beamlet_number: str
+ :param input_data: input data to be processed
+ :type input_data: dict(dict(numpy.ndarray))
+ :return:
+ """
+ normalized_data = dict()
+ for station in input_data:
+ normalized_data_per_station = dict()
+ normalized_data[station] = normalized_data_per_station
+
+ data_per_station = input_data[station]
+
+ for frequency in data_per_station:
+ data_per_station_per_frequency = data_per_station[frequency]
+
+ normalized_data_per_station_frequency = \
+ normalize_crosscorrelation_beams(data_per_station_per_frequency,
+ central_beamlet_number)
+
+ normalized_data_per_station[frequency] = normalized_data_per_station_frequency
+ return normalized_data
+
+
+def normalize_crosscorrelation_beams(data_per_station_frequency,
+ central_beamlet):
+ """
+ It normalize the crosscorrelation per each beam given the normalization array computed
+ from the central beam
+ :param data_per_station_frequency:
+ :param normalization_array:
+ :return:
+ """
+
+ """
+ TO avoid copying the data in a different array and then copy the data back
+ I compute directly the product of the two array of matrices
+
+ the resulting cross polarization as a function of time are:
+ XX = XX_b * XX_n + XY_b * YX_n
+ XY = XX_b * XY_n + XY_b * YY_n
+ YX = YX_b * XX_n + YY_b * YX_n
+ YY = YX_b * XY_n + YY_b * YY_n
+
+ where b is the beam and n is the normalization array
+ """
+ is_invertion_matrix_valid = True
+ normalization_array = None
+ try:
+ normalization_array = invert_central_beam_crosscorrelations(data_per_station_frequency,
+ central_beamlet)
+ except numpy.linalg.LinAlgError as e:
+ logger.warning('error inverting central beamlet crosscorrelation for beamlet %s: %s',
+ central_beamlet,
+ e)
+ print('error inverting central beamlet crosscorrelation for beamlet %s: %s' % (
+ central_beamlet,
+ e))
+ is_invertion_matrix_valid = False
+
+ normalized_data_per_beam = dict()
+ for beam in data_per_station_frequency:
+
+ data_per_beam = data_per_station_frequency[beam]
+ normalized_data = numpy.array(data_per_beam)
+ normalized_data_per_beam[beam] = normalized_data
+
+ if is_invertion_matrix_valid == False:
+ normalized_data['flag'] = True
+ continue
+
+ normalized_data['XX'] = \
+ data_per_beam['XX'] * normalization_array[__POLARIZATION_TO_INDEX['XX']] + \
+ data_per_beam['XY'] * normalization_array[__POLARIZATION_TO_INDEX['YX']]
+ normalized_data['XY'] = \
+ data_per_beam['XX'] * normalization_array[__POLARIZATION_TO_INDEX['XY']] + \
+ data_per_beam['XY'] * normalization_array[__POLARIZATION_TO_INDEX['YY']]
+ normalized_data['YX'] = \
+ data_per_beam['YX'] * normalization_array[__POLARIZATION_TO_INDEX['XX']] + \
+ data_per_beam['YY'] * normalization_array[__POLARIZATION_TO_INDEX['YX']]
+ normalized_data['YY'] = \
+ data_per_beam['YX'] * normalization_array[__POLARIZATION_TO_INDEX['XY']] + \
+ data_per_beam['YY'] * normalization_array[__POLARIZATION_TO_INDEX['YY']]
+ return normalized_data_per_beam
+
+
+def invert_central_beam_crosscorrelations(data_per_station_per_frequency, central_beamlet_number):
+ data_per_central_beamlet = __extract_crosscorrelation_matrices_from_data(
+ data_per_station_per_frequency[central_beamlet_number])
+
+ return __invert_crosscorrelation_matrices(data_per_central_beamlet)
+
+
+def __extract_crosscorrelation_matrices_from_data(data_per_station_frequency_beam):
+ new_shape = (data_per_station_frequency_beam.shape[0], 2, 2)
+ new_array = numpy.zeros(shape=new_shape, dtype=numpy.complex)
+
+ new_array[:, 0, 0] = data_per_station_frequency_beam['XX']
+ new_array[:, 0, 1] = data_per_station_frequency_beam['XY']
+ new_array[:, 1, 0] = data_per_station_frequency_beam['YX']
+ new_array[:, 1, 1] = data_per_station_frequency_beam['YY']
+
+ return new_array
+
+
+def __invert_crosscorrelation_matrices(cross_correlation_matrices):
+ """
+ Invert the matrices along the last axis
+ :param cross_correlation_matrices: matrices to invert
+ :type cross_correlation_matrices: numpy.ndarray
+ :return:
+ """
+ assert cross_correlation_matrices.ndim >= 3
+
+ result_array = numpy.rollaxis(numpy.linalg.inv(cross_correlation_matrices), 0, 3)
+
+ return result_array
+
diff --git a/CAL/CalibrationProcessing/lib/processing/solver.py b/CAL/CalibrationProcessing/lib/processing/solver.py
new file mode 100644
index 0000000000000000000000000000000000000000..9068eda32d97c0010effd695fb83a1b93a94b195
--- /dev/null
+++ b/CAL/CalibrationProcessing/lib/processing/solver.py
@@ -0,0 +1,169 @@
+import logging
+
+import lofar.calibration.common.datacontainers as datacontainers
+from scipy.constants import c as LIGHT_SPEED
+
+logger = logging.getLogger(__name__)
+
+import numpy
+
+__FRQ_TO_PHASE = 2. * numpy.pi / LIGHT_SPEED
+
+
+def _rotate_antenna_coordinates(dataset):
+ """
+
+ :param dataset:
+ :type dataset: datacontainers.HolographyDataset
+ :return:
+ """
+ station_position = numpy.array(dataset.target_station_position)
+ antenna_field_position = numpy.array(dataset.antenna_field_position)
+
+ antenna_position_offset = numpy.array([station_position - antenna_position for
+ antenna_position in antenna_field_position])
+
+ station_rotation_matrix = dataset.rotation_matrix
+ return numpy.dot(station_rotation_matrix, antenna_position_offset.T).T
+
+
+def _compute_expected_phase_delay(l, m, x, y, frequency):
+ """
+ Convert the antenna position (x,y) and the pointing cosines (l,m) into a phase delay for
+ the given frequency
+ :param l:
+ :param m:
+ :param x:
+ :param y:
+ :param frequency:
+ :return:
+ """
+ phase = (x * l + y * m) * __FRQ_TO_PHASE * frequency
+ return phase
+
+
+def _compute_pointing_matrices_per_station_frequency_beam(dataset, datatable, frequency):
+ """
+
+ :param dataset:
+ :type dataset: datacontainers.HolographyDataset
+ :param datatable:
+ :type table to reduce:
+ :param frequency:
+ :type frequency: float
+ :return:
+ """
+
+ rotated_coordinates = _rotate_antenna_coordinates(dataset)
+ n_antennas = len(dataset.antenna_field_position)
+ one_over_n_antennas = 1. / float(n_antennas)
+ n_beams = len(datatable)
+
+ pointing_matrix = numpy.matrix(numpy.zeros((n_beams, n_antennas), dtype=complex))
+
+ for i in range(n_beams):
+ for j in range(n_antennas):
+ l = datatable[str(i)]['mean']['l'][0]
+ m = datatable[str(i)]['mean']['m'][0]
+ x, y = rotated_coordinates[j, 0:2]
+ phase = _compute_expected_phase_delay(l, m, x, y, frequency)
+ pointing_matrix[i, j] = numpy.exp(-1.j * phase) * one_over_n_antennas
+
+ return pointing_matrix
+
+
+def _solve_gains_complex(visibilities, matrix):
+ """
+ To solve a complex linear system of equation it is possible to rewrite it in a equivalent
+ real system.
+
+ [re1 + i + im1 , re2 + i + im2 ; [ re1, -im1, re2, -im2]
+ ~ [ im1, re1, im2, re2]
+ re3 + i + im3 , re4 + i + im4 ] [ re3, -im3, re4, -im4]
+ [ re3, re3, re$, -im4]
+
+ :param matrix:
+ :param visibilities:
+ :return:
+ """
+ matrix_real = numpy.matrix(numpy.zeros((matrix.shape[0]*2, matrix.shape[1]*2)))
+ matrix_real[0::2, 0::2] = matrix.real
+ matrix_real[0::2, 1::2] = -matrix.imag
+ matrix_real[1::2, 0::2] = matrix.real
+ matrix_real[1::2, 1::2] = matrix.imag
+
+ visibilities_real = numpy.matrix(numpy.zeros(visibilities.shape[0] * 2)).T
+ visibilities_real[0::2] = visibilities.real
+ visibilities_real[1::2] = visibilities.imag
+
+ result = _solve_gains(visibilities_real, matrix_real)
+ if result['flag'] == False:
+ noise = result['relative_error']
+
+ noise = numpy.sqrt(noise[0::2].A1**2. + noise[1::2].A1**2.)
+ result["relative_error"] = noise
+ return result
+
+
+def _solve_gains(visibilities, matrix):
+
+ try:
+ gains, residual, _, _ = numpy.linalg.lstsq(matrix, visibilities, rcond=None)
+ noise = abs(matrix * gains - visibilities) / abs(visibilities)
+
+ return dict(gains=gains, residual=residual, relative_error=noise,
+ flag=False)
+ except Exception as e:
+ logger.warning('error solving for the gains: %s', e)
+ __empty = dict(gains=numpy.array(numpy.nan),
+ residual=numpy.array(numpy.nan),
+ relative_error=numpy.array(numpy.nan),
+ flag=numpy.array(True))
+ return __empty
+
+
+def _solve_gains_least_square(dataset, datatable, frequency, direct_complex=False):
+ """
+ SOLVE THE EQUATION M * G = V FOR G
+ :param dataset:
+ :type dataset: datacontainers.HolographyDataset
+ :param datatable:
+ :param frequency:
+ :type frequency: float
+ :return:
+ """
+ matrix = _compute_pointing_matrices_per_station_frequency_beam(dataset, datatable, frequency)
+ n_beams = len(datatable)
+
+ result = dict()
+ for polarization in ['XX', 'XY', 'YX', 'YY']:
+ visibilities = numpy.matrix([datatable[str(i)]['mean'][polarization] for i in range(n_beams)])
+ if direct_complex == True:
+ result[polarization] = _solve_gains(visibilities, matrix)
+ else:
+ result[polarization] = _solve_gains_complex(visibilities, matrix)
+
+ return result
+
+
+def solve_gains_per_datatable(dataset, datatable, **kwargs):
+ """
+
+ :param dataset:
+ :type dataset: datacontainers.HolographyDataset
+ :param datatable:
+ :type dataset: dict(dict(numpy.ndarray))
+ :return:
+ """
+ result = dict()
+
+ for station in datatable:
+ result_per_station = dict()
+ result[station] = result_per_station
+ data_per_station = datatable[station]
+ for frequency in data_per_station:
+ data_per_frequency = data_per_station[frequency]
+ result_per_station[str(frequency)] = \
+ _solve_gains_least_square(dataset, data_per_frequency, float(frequency), **kwargs)
+
+ return result
diff --git a/CAL/CalibrationProcessing/test/CMakeLists.txt b/CAL/CalibrationProcessing/test/CMakeLists.txt
index e9e6ce8708ff4c0823725d09dc7040b96fd3dac0..7a639e742fde607a7c5f9c397986e1acec82bbd8 100644
--- a/CAL/CalibrationProcessing/test/CMakeLists.txt
+++ b/CAL/CalibrationProcessing/test/CMakeLists.txt
@@ -19,6 +19,11 @@
include(LofarCTest)
-lofar_add_test(t_processing)
+python_install(
+ __init__.py
+ utils.py
+ DESTINATION lofar/calibration/testing)
+
+add_subdirectory(processing)
diff --git a/CAL/CalibrationProcessing/test/__init__.py b/CAL/CalibrationProcessing/test/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/CAL/CalibrationProcessing/test/processing/CMakeLists.txt b/CAL/CalibrationProcessing/test/processing/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e9e6ce8708ff4c0823725d09dc7040b96fd3dac0
--- /dev/null
+++ b/CAL/CalibrationProcessing/test/processing/CMakeLists.txt
@@ -0,0 +1,24 @@
+# Copyright (C) 2018 ASTRON (Netherlands Institute for Radio Astronomy)
+# P.O. Box 2, 7990 AA Dwingeloo, The Netherlands
+#
+# This file is part of the LOFAR software suite.
+# The LOFAR software suite is free software: you can redistribute it and/or
+# modify it under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# The LOFAR software suite is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
+
+# $Id$
+include(LofarCTest)
+
+
+lofar_add_test(t_processing)
+
+
diff --git a/CAL/CalibrationProcessing/test/t_processing.py b/CAL/CalibrationProcessing/test/processing/t_processing.py
similarity index 64%
rename from CAL/CalibrationProcessing/test/t_processing.py
rename to CAL/CalibrationProcessing/test/processing/t_processing.py
index 83c2f1bbd1e78dc529d626dd73337e99ee113dc7..3319852cf0c502e5971f3dd8c02905ae8fc3a973 100755
--- a/CAL/CalibrationProcessing/test/t_processing.py
+++ b/CAL/CalibrationProcessing/test/processing/t_processing.py
@@ -1,21 +1,72 @@
import logging
import unittest
+import os
from lofar.calibration.common.datacontainers import HolographyDataset
from lofar.calibration.processing import average_dataset_by_sample, average_values_by_sample
from lofar.calibration.common.datacontainers.holography_dataset_definitions import HDS_data_sample_type
-from lofar.calibration.processing.matrices import extract_crosscorrelation_matrices_from_HDS_datatable
+from lofar.calibration.processing.normalize import __extract_crosscorrelation_matrices_from_data
+
+from lofar.calibration.testing.utils import create_holography_dataset
+import lofar.calibration.testing.utils as utils
logger = logging.getLogger('t_processing')
import numpy
+import matplotlib.pyplot as plt
# READ doc/Holography_Data_Set.md! It contains the location from which the
# test data must be downloaded and this file generated.
SAMPLE_DATASET = '/data/test/HolographyObservation/CS001HBA0.hdf5'
+def is_sample_dataset_available():
+ return os.path.exists(SAMPLE_DATASET)
+
+
+class TestSignalGenerator(unittest.TestCase):
+ def test_gaussian_noise_generation(self):
+ test_average_value = 5. + 12j
+ test_real_std = 12
+ test_imag_std = 3
+ test_sample_size = 100
+ gaussian_generator = utils.TestComplexGaussianNoise(test_average_value,
+ test_real_std,
+ test_imag_std)
+ sample_noise = gaussian_generator.generate(test_sample_size, seed=5)
+ self.assertGreater(numpy.std(sample_noise.imag)/test_imag_std, .9)
+ self.assertGreater(numpy.std(sample_noise.imag)/test_imag_std, .9)
+ self.assertGreater(abs(numpy.average(sample_noise).real / test_average_value.real), .9)
+ self.assertGreater(abs(numpy.average(sample_noise).imag / test_average_value.imag), .9)
+
+ def test_signal_generator(self):
+ test_imag_coeff = 3.j
+ test_real_coeff = 2.
+ expected_dependencies_imag_real = test_imag_coeff/test_real_coeff
+ linear_function = lambda x: 2 * x + 3.j * x
+
+ test_average_value = 5. + 12j
+ test_real_std = 12
+ test_imag_std = 3
+ test_time_samples = numpy.linspace(0, 100, 1000).tolist()
+ gaussian_generator = utils.TestComplexGaussianNoise(test_average_value,
+ test_real_std,
+ test_imag_std)
+
+ signal_generator = utils.TestSignal(linear_function, gaussian_generator)
+ sample_signal = signal_generator.generate(test_time_samples)
+
+ coefficients = numpy.polyfit(sample_signal.real, sample_signal.imag, 1)
+ m, q = coefficients[0], coefficients[1]
+ self.assertGreater(abs(m/expected_dependencies_imag_real), .9 )
+
class TestProcessing(unittest.TestCase):
+ sample_dataset = None
def setUp(self):
- pass
+ self.sample_dataset = create_holography_dataset()
+
+ def test_sample_dataset(self):
+ for reference_station in self.sample_dataset.reference_stations:
+ self.assertIn(reference_station, self.sample_dataset.data)
+
def test_average_array_by_sample_real_values(self):
"""
@@ -51,6 +102,7 @@ class TestProcessing(unittest.TestCase):
numpy.testing.assert_almost_equal(result['std_real'], expected_result_std.real)
numpy.testing.assert_almost_equal(result['std_imag'], expected_result_std.real)
+ @unittest.skipIf(is_sample_dataset_available(), 'missing test dataset')
def test_average_datatable_by_sample(self):
holography_dataset = HolographyDataset.load_from_file(SAMPLE_DATASET)
averaged_data = average_dataset_by_sample(SAMPLE_DATASET, 2).derived_data[
@@ -85,7 +137,7 @@ class TestProcessing(unittest.TestCase):
test_data['YX'] = numpy.random.randn(3) + 1.j * numpy.random.randn(3)
test_data['YY'] = numpy.random.randn(3) + 1.j * numpy.random.randn(3)
- matrices_from_data = extract_crosscorrelation_matrices_from_HDS_datatable(test_data)
+ matrices_from_data = __extract_crosscorrelation_matrices_from_data(test_data)
numpy.testing.assert_almost_equal(test_data['XX'], matrices_from_data[0, 0, :])
numpy.testing.assert_almost_equal(test_data['XY'], matrices_from_data[0, 1, :])
diff --git a/CAL/CalibrationProcessing/test/t_processing.run b/CAL/CalibrationProcessing/test/processing/t_processing.run
similarity index 100%
rename from CAL/CalibrationProcessing/test/t_processing.run
rename to CAL/CalibrationProcessing/test/processing/t_processing.run
diff --git a/CAL/CalibrationProcessing/test/t_processing.sh b/CAL/CalibrationProcessing/test/processing/t_processing.sh
similarity index 100%
rename from CAL/CalibrationProcessing/test/t_processing.sh
rename to CAL/CalibrationProcessing/test/processing/t_processing.sh
diff --git a/CAL/CalibrationProcessing/test/profiling_processing.py b/CAL/CalibrationProcessing/test/profiling_processing.py
index 0ecbcf3989c3b5a40bb671b496b457b3b3f603ce..ea0bbaf9fa01bb181c9f7d0f69db94f712684b3e 100644
--- a/CAL/CalibrationProcessing/test/profiling_processing.py
+++ b/CAL/CalibrationProcessing/test/profiling_processing.py
@@ -54,6 +54,7 @@ class profile_averaging_per_sample(ProfilingTest):
def process(self):
processing.average_dataset_by_sample(SAMPLE_DATASET, 2)
+
if __name__ == '__main__':
profile_holography_dataset_read().run().print_results()
profile_averaging_per_sample().run().print_results()
diff --git a/CAL/CalibrationProcessing/test/test_utils.py b/CAL/CalibrationProcessing/test/test_utils.py
deleted file mode 100644
index ee96891cef17a8a7473da52f59de6df3b7bdafc6..0000000000000000000000000000000000000000
--- a/CAL/CalibrationProcessing/test/test_utils.py
+++ /dev/null
@@ -1,16 +0,0 @@
-
-from lofar.calibration.common.datacontainers.holography_dataset import HolographyDataset
-
-
-def create_holography_dataset(target_station='CS001',
- reference_stations=('RS301', 'RS302'),
- pointings=(
-
- )
- ):
- """
- Creates a mock Holography Dataset that can be used for testing
- :param target_station:
- :param reference_stations:
- :return:
- """
\ No newline at end of file
diff --git a/CAL/CalibrationProcessing/test/utils.py b/CAL/CalibrationProcessing/test/utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..d0696d0fde2a9221ff18a0126626cbbc42e7fd9c
--- /dev/null
+++ b/CAL/CalibrationProcessing/test/utils.py
@@ -0,0 +1,255 @@
+import numpy
+from collections import defaultdict
+
+from lofar.calibration.common.datacontainers.holography_dataset import HolographyDataset
+from lofar.calibration.common.datacontainers.holography_dataset_definitions import HDS_coordinate_type, HDS_data_sample_type
+
+DEFAULT_TARGET_STATION = 'CS001'
+DEFAULT_TARGET_STATION_POSITION = numpy.array([0, 0, 0])
+
+DEFAULT_REFERENCE_STATIONS = ('RS301', 'RS302')
+DEFAULT_SOURCE_NAME = 'test source'
+DEFAULT_SOURCE_POINTING = numpy.array([(1.5, 0.9, 'J2000')], HDS_coordinate_type)
+DEFAULT_SAS_IDS = numpy.array([123456])
+DEFAULT_SOURCE_POSITION = numpy.array([3826896.19174 , 460979.502113, 5064658.231])
+DEFAULT_BEAMLETS = numpy.array([0, 1, 2, 3 , 4])
+DEFAULT_FREQUENCIES = [150000.0]
+DEFAULT_POINTINGS = {str(DEFAULT_FREQUENCIES[0]): {
+ '0': numpy.array([(1.5, 0.9, 'J2000')], HDS_coordinate_type),
+ '1': numpy.array([(2.0, 0.9, 'J2000')], HDS_coordinate_type),
+ '2': numpy.array([(1.0, 0.9, 'J2000')], HDS_coordinate_type),
+ '3': numpy.array([(1.0, 0.8, 'J2000')], HDS_coordinate_type),
+ '4': numpy.array([(1.0, 1, 'J2000')], HDS_coordinate_type)}}
+
+DEFAULT_ANTENNA_FIELD_POSITION = numpy.array([[1, 0, 1],
+ [-1, 0, 1],
+ [+1, 0, 1]])
+
+DEFAULT_ROTATION_MATRIX = numpy.array([[-1.19595e-01, 9.92823e-01, 3.30000e-05],
+ [-7.91954e-01, -9.54190e-02, 6.03078e-01],
+ [ 5.98753e-01, 7.20990e-02, 7.97682e-01]])
+
+DEFAULT_MODE = 5
+DEFAULT_TIME_SAMPLE = [1, 2, 3, 4, 5, 6]
+
+
+class TestComplexGaussianNoise():
+ def __init__(self, average, std_real, std_imag):
+ """
+ :param average: average value
+ :type average: complex
+ :param std_real: standard deviation for the real part
+ :type std_real: float
+ :param std_imag: standard deviation for the imaginary part
+ :type std_imag: float
+ """
+ self.average = average
+ self.std_real = std_real
+ self.std_imag = std_imag
+
+ def generate(self, sample_size, seed=None):
+ """
+ Generate a series of random values
+ :param sample_size:
+ :type sample_size: int
+ :param seed:
+ :type seed: int
+ :return:
+ """
+ if seed:
+ numpy.random.seed(seed)
+
+ real = self.std_real * numpy.random.randn(sample_size)
+ imag = self.std_imag * numpy.random.randn(sample_size) * 1.j
+
+ return real + imag + self.average
+
+
+class TestSignal():
+ def __init__(self, function, noise):
+ """
+
+ :param function:
+ :type function: function
+ :param noise:
+ :type noise: TestComplexGaussianNoise
+ """
+ self.function = function
+ self.noise = noise
+
+ def generate(self, timesamples, seed=None):
+ """
+
+ :param timesamples: list of time samples in which to compute the signal
+ :type timesamples: Iterable[float]
+ :param seed: seed for the gaussian signal generator
+ :type seed: None|int
+ :return: the generated sample signal
+ """
+ n_samples = len(timesamples)
+ result = numpy.zeros(n_samples, dtype=complex)
+ for i, timesample in enumerate(timesamples):
+ result[i] = self.function(timesample)
+ result += self.noise.generate(n_samples, seed)
+
+ return result
+
+
+DEFAULT_AVERAGE_VALUES_PER_STATION_AND_POINTING = {
+ "RS301": {
+ '150000.0': {
+ '0': {
+ "mean": 1+5j,
+ "std_imag": 1,
+ "std_real": 2,
+ "real_coeff": 0,
+ "imag_coeff": 0
+ },
+ '1': {
+ "mean": 1 + 5j,
+ "std_imag": 1,
+ "std_real": 2,
+ "real_coeff": 0,
+ "imag_coeff": 0
+ },
+ '2': {
+ "mean": 1 + 5j,
+ "std_imag": 1,
+ "std_real": 2,
+ "real_coeff": 0,
+ "imag_coeff": 0
+ },
+ '3': {
+ "mean": 1 + 5j,
+ "std_imag": 1,
+ "std_real": 2,
+ "real_coeff": 0,
+ "imag_coeff": 0
+ },
+ '4': {
+ "mean": 1 + 5j,
+ "std_imag": 1,
+ "std_real": 2,
+ "real_coeff": 0,
+ "imag_coeff": 0
+ }
+ }
+ },
+ "RS302": {
+ '150000.0': {
+ '0': {
+ "mean": 1+5j,
+ "std_imag": 1,
+ "std_real": 2,
+ "real_coeff": 0,
+ "imag_coeff": 0
+ },
+ '1': {
+ "mean": 1 + 5j,
+ "std_imag": 1,
+ "std_real": 2,
+ "real_coeff": 0,
+ "imag_coeff": 0
+ },
+ '2': {
+ "mean": 1 + 5j,
+ "std_imag": 1,
+ "std_real": 2,
+ "real_coeff": 0,
+ "imag_coeff": 0
+ },
+ '3': {
+ "mean": 1 + 5j,
+ "std_imag": 1,
+ "std_real": 2,
+ "real_coeff": 0,
+ "imag_coeff": 0
+ },
+ '4': {
+ "mean": 1 + 5j,
+ "std_imag": 1,
+ "std_real": 2,
+ "real_coeff": 0,
+ "imag_coeff": 0
+ }
+ }
+ }
+}
+
+
+def recursive_default_dict_generator():
+ return defaultdict(recursive_default_dict_generator)
+
+
+def generate_single_pointing_dataset(timesamples, parameters):
+ """
+
+ :param timesamples:
+ :type timesamples: list
+ :return:
+ """
+ dataset = numpy.zeros(len(timesamples), dtype=HDS_data_sample_type)
+ return dataset
+
+
+def generate_sample_data(holography_dataset,
+ timesamples = DEFAULT_TIME_SAMPLE,
+ data_statistical_characteristics=DEFAULT_AVERAGE_VALUES_PER_STATION_AND_POINTING):
+ """
+ Generates a sample data to start testing the processing step
+ :param holography_dataset: dataset prefilled with the header parameters
+ :type holography_dataset: HolographyDataset
+ :param timesamples:
+
+ :return:
+ """
+ holography_dataset.data = recursive_default_dict_generator()
+
+ for reference_station in holography_dataset.reference_stations:
+ for frequency in holography_dataset.frequencies:
+ frequency_string = str(frequency)
+ for beamlet in holography_dataset.beamlets:
+ beamlet_string = str(beamlet)
+ holography_dataset.data[reference_station][frequency_string][beamlet_string] =\
+ generate_single_pointing_dataset(timesamples,
+ data_statistical_characteristics[reference_station][frequency_string][beamlet_string])
+
+
+
+
+def create_holography_dataset(target_station_name=DEFAULT_TARGET_STATION,
+ target_station_position = DEFAULT_TARGET_STATION_POSITION,
+ reference_stations=DEFAULT_REFERENCE_STATIONS,
+ timesamples = DEFAULT_TIME_SAMPLE,
+ pointings=DEFAULT_POINTINGS,
+ sas_ids = DEFAULT_SAS_IDS,
+ source_name=DEFAULT_SOURCE_NAME,
+ source_position=DEFAULT_SOURCE_POINTING,
+ mode=DEFAULT_MODE,
+ beamlets=DEFAULT_BEAMLETS,
+ antenna_field_positions=DEFAULT_ANTENNA_FIELD_POSITION,
+ frequencies=DEFAULT_FREQUENCIES,
+ data_statistical_characteristics=DEFAULT_AVERAGE_VALUES_PER_STATION_AND_POINTING
+ ):
+ """
+ Creates a mock Holography Dataset that can be used for testing
+ """
+ dataset = HolographyDataset()
+ dataset.mode = mode
+ dataset.target_station_position = target_station_position
+ dataset.target_station_name = target_station_name
+ dataset.beamlets = beamlets
+ dataset.source_position = source_position
+ dataset.reference_stations = reference_stations
+ dataset.antenna_field_position = antenna_field_positions
+ dataset.source_name = source_name
+ dataset.sas_ids = sas_ids
+ dataset.ra_dec = pointings
+ dataset.central_beamlets = dataset.find_central_beamlets(source_position,
+ pointings,
+ frequencies,
+ beamlets)
+ dataset.frequencies = frequencies
+
+ generate_sample_data(dataset, timesamples)
+ return dataset