diff --git a/examples/caltable_example.ipynb b/examples/caltable_example.ipynb index 4d8c0db..bf94734 100644 --- a/examples/caltable_example.ipynb +++ b/examples/caltable_example.ipynb @@ -2,10 +2,13 @@ "cells": [ { "cell_type": "code", - "execution_count": 2, "id": "6e4714df-ebca-4a58-bb19-c2e931f168f4", - "metadata": {}, - "outputs": [], + "metadata": { + "ExecuteTime": { + "end_time": "2025-02-04T15:37:54.839002Z", + "start_time": "2025-02-04T15:37:51.528062Z" + } + }, "source": [ "import toolviper\n", "import calviper\n", @@ -17,25 +20,51 @@ "\n", "from xradio.measurement_set.convert_msv2_to_processing_set import convert_msv2_to_processing_set\n", "from xradio.measurement_set import open_processing_set" - ] + ], + "outputs": [], + "execution_count": 1 }, { "cell_type": "code", - "execution_count": 4, "id": "24a3971d-c7d5-4aef-92dc-d90bf2750498", - "metadata": {}, + "metadata": { + "ExecuteTime": { + "end_time": "2025-02-04T15:39:25.015619Z", + "start_time": "2025-02-04T15:39:22.082556Z" + } + }, + "source": [ + "toolviper.utils.data.download(file=\"Antennae_North.cal.lsrk.split.ms\")\n", + "\n", + "_ = convert_msv2_to_processing_set(\n", + " in_file=\"Antennae_North.cal.lsrk.split.ms\",\n", + " out_file=\"Antennae_North.cal.lsrk.split.vis.zarr\",\n", + " parallel=False,\n", + " overwrite=True,\n", + " main_chunksize={\n", + " \"frequency\": 3\n", + " },\n", + ")" + ], "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "[\u001B[38;2;128;05;128m2024-11-21 09:29:37,320\u001B[0m] \u001B[38;2;255;160;0m WARNING\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m File exists: \u001B[38;2;50;50;205m/export/home/ajax/jhoskins/.conda/envs/calviper/lib/python3.10/site-packages/toolviper/utils/data/.dropbox\u001B[0m \n", - "[\u001B[38;2;128;05;128m2024-11-21 09:29:37,321\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m Updating file metadata information ... \n", + "[\u001B[38;2;128;05;128m2025-02-04 10:39:22,085\u001B[0m] \u001B[38;2;255;160;0m WARNING\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m File exists: \u001B[38;2;50;50;205m/home/mystletainn/miniconda3/envs/calviper/lib/python3.10/site-packages/toolviper/utils/data/.dropbox\u001B[0m \n", + "[\u001B[38;2;128;05;128m2025-02-04 10:39:22,086\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m Updating file metadata information ... \n", " " ] }, { "data": { + "text/plain": [ + " \n", + " \u001B[1m \u001B[0m\u001B[1mDownload List \u001B[0m\u001B[1m \u001B[0m \n", + " ────────────────────────────────── \n", + " \u001B[35mAntennae_North.cal.lsrk.split.ms\u001B[0m \n", + " \n" + ], "text/html": [ "
\n", " Download List \n", @@ -43,13 +72,6 @@ " Antennae_North.cal.lsrk.split.ms \n", " \n", "\n" - ], - "text/plain": [ - " \n", - " \u001B[1m \u001B[0m\u001B[1mDownload List \u001B[0m\u001B[1m \u001B[0m \n", - " ────────────────────────────────── \n", - " \u001B[35mAntennae_North.cal.lsrk.split.ms\u001B[0m \n", - " \n" ] }, "metadata": {}, @@ -59,46 +81,36 @@ "name": "stdout", "output_type": "stream", "text": [ - "[\u001B[38;2;128;05;128m2024-11-21 09:29:39,031\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m File exists: Antennae_North.cal.lsrk.split.ms \n", - "[\u001B[38;2;128;05;128m2024-11-21 09:29:39,053\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m Partition scheme that will be used: ['DATA_DESC_ID', 'OBS_MODE', 'OBSERVATION_ID', 'FIELD_ID'] \n", - "[\u001B[38;2;128;05;128m2024-11-21 09:29:39,248\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m Number of partitions: 12 \n", - "[\u001B[38;2;128;05;128m2024-11-21 09:29:39,250\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [0], DDI [0], STATE [32 23 30 37], FIELD [0], SCAN [ 9 17 21 25] \n", - "[\u001B[38;2;128;05;128m2024-11-21 09:29:39,474\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [0], DDI [0], STATE [33 24 31], FIELD [1], SCAN [ 9 17 21] \n", - "[\u001B[38;2;128;05;128m2024-11-21 09:29:39,625\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [0], DDI [0], STATE [34 25 32], FIELD [2], SCAN [ 9 17 21] \n", - "[\u001B[38;2;128;05;128m2024-11-21 09:29:39,772\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [1], DDI [0], STATE [32 23 30 37], FIELD [0], SCAN [26 34 38 42] \n", - "[\u001B[38;2;128;05;128m2024-11-21 09:29:39,940\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [1], DDI [0], STATE [33 24 31], FIELD [1], SCAN [26 34 38] \n", - "[\u001B[38;2;128;05;128m2024-11-21 09:29:40,107\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [1], DDI [0], STATE [34 25 32], FIELD [2], SCAN [26 34 38] \n", - "[\u001B[38;2;128;05;128m2024-11-21 09:29:40,265\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [2], DDI [0], STATE [32], FIELD [0], SCAN [43] \n", - "[\u001B[38;2;128;05;128m2024-11-21 09:29:40,452\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [2], DDI [0], STATE [33], FIELD [1], SCAN [43] \n", - "[\u001B[38;2;128;05;128m2024-11-21 09:29:40,624\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [2], DDI [0], STATE [34], FIELD [2], SCAN [43] \n", - "[\u001B[38;2;128;05;128m2024-11-21 09:29:40,790\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [3], DDI [0], STATE [48 39 46 53], FIELD [0], SCAN [48 56 60 64] \n", - "[\u001B[38;2;128;05;128m2024-11-21 09:29:40,953\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [3], DDI [0], STATE [49 40 47], FIELD [1], SCAN [48 56 60] \n", - "[\u001B[38;2;128;05;128m2024-11-21 09:29:41,114\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [3], DDI [0], STATE [50 41 48], FIELD [2], SCAN [48 56 60] \n" + "[\u001B[38;2;128;05;128m2025-02-04 10:39:22,790\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m File exists: Antennae_North.cal.lsrk.split.ms \n", + "[\u001B[38;2;128;05;128m2025-02-04 10:39:22,842\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m Partition scheme that will be used: ['DATA_DESC_ID', 'OBS_MODE', 'OBSERVATION_ID', 'FIELD_ID'] \n", + "[\u001B[38;2;128;05;128m2025-02-04 10:39:22,932\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m Number of partitions: 12 \n", + "[\u001B[38;2;128;05;128m2025-02-04 10:39:22,933\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [0], DDI [0], STATE [32 23 30 37], FIELD [0], SCAN [ 9 17 21 25] \n", + "[\u001B[38;2;128;05;128m2025-02-04 10:39:23,108\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [0], DDI [0], STATE [33 24 31], FIELD [1], SCAN [ 9 17 21] \n", + "[\u001B[38;2;128;05;128m2025-02-04 10:39:23,244\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [0], DDI [0], STATE [34 25 32], FIELD [2], SCAN [ 9 17 21] \n", + "[\u001B[38;2;128;05;128m2025-02-04 10:39:23,379\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [1], DDI [0], STATE [32 23 30 37], FIELD [0], SCAN [26 34 38 42] \n", + "[\u001B[38;2;128;05;128m2025-02-04 10:39:23,548\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [1], DDI [0], STATE [33 24 31], FIELD [1], SCAN [26 34 38] \n", + "[\u001B[38;2;128;05;128m2025-02-04 10:39:23,726\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [1], DDI [0], STATE [34 25 32], FIELD [2], SCAN [26 34 38] \n", + "[\u001B[38;2;128;05;128m2025-02-04 10:39:23,909\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [2], DDI [0], STATE [32], FIELD [0], SCAN [43] \n", + "[\u001B[38;2;128;05;128m2025-02-04 10:39:24,091\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [2], DDI [0], STATE [33], FIELD [1], SCAN [43] \n", + "[\u001B[38;2;128;05;128m2025-02-04 10:39:24,259\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [2], DDI [0], STATE [34], FIELD [2], SCAN [43] \n", + "[\u001B[38;2;128;05;128m2025-02-04 10:39:24,449\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [3], DDI [0], STATE [48 39 46 53], FIELD [0], SCAN [48 56 60 64] \n", + "[\u001B[38;2;128;05;128m2025-02-04 10:39:24,623\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [3], DDI [0], STATE [49 40 47], FIELD [1], SCAN [48 56 60] \n", + "[\u001B[38;2;128;05;128m2025-02-04 10:39:24,840\u001B[0m] \u001B[38;2;50;50;205m INFO\u001B[0m\u001B[38;2;112;128;144m viperlog: \u001B[0m OBSERVATION_ID [3], DDI [0], STATE [50 41 48], FIELD [2], SCAN [48 56 60] \n" ] } ], - "source": [ - "toolviper.utils.data.download(file=\"../Antennae_North.cal.lsrk.split.ms\")\n", - "\n", - "_ = convert_msv2_to_processing_set(\n", - " in_file=\"../Antennae_North.cal.lsrk.split.ms\",\n", - " out_file=\"../Antennae_North.cal.lsrk.split.vis.zarr\",\n", - " parallel=False,\n", - " overwrite=True,\n", - " main_chunksize={\n", - " \"frequency\": 3\n", - " },\n", - ")" - ] + "execution_count": 5 }, { "cell_type": "code", - "execution_count": 6, "id": "4077e188-728a-4c90-9d31-812b9c888756", "metadata": { - "scrolled": true + "scrolled": true, + "ExecuteTime": { + "end_time": "2025-02-04T15:39:27.159422Z", + "start_time": "2025-02-04T15:39:26.828846Z" + } }, - "outputs": [], "source": [ "ps = open_processing_set(\n", " ps_store=\"Antennae_North.cal.lsrk.split.vis.zarr\",\n", @@ -106,42 +118,48 @@ ")\n", "\n", "#ps.summary()" - ] + ], + "outputs": [], + "execution_count": 6 }, { "cell_type": "code", - "execution_count": 7, "id": "eebfc16c-13ec-4a31-830c-5cac52d77863", - "metadata": {}, - "outputs": [], - "source": [ - "mxds = ps[\"Antennae_North.cal.lsrk.split_00\"]" - ] + "metadata": { + "ExecuteTime": { + "end_time": "2025-02-04T17:08:45.691150Z", + "start_time": "2025-02-04T17:08:45.684895Z" + } + }, + "source": "mxds = ps[\"Antennae_North.cal.lsrk.split_00\"]", + "outputs": [ + { + "data": { + "text/plain": [ + "['X', 'X', 'Y', 'Y']" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "execution_count": 12 }, { "cell_type": "code", - "execution_count": 8, "id": "166e3370-c008-44e0-915d-2e1304c228df", - "metadata": {}, - "outputs": [ - { - "ename": "TypeError", - "evalue": "CalibrationTable() missing 1 required positional argument: 'self'", - "output_type": "error", - "traceback": [ - "\u001B[0;31m---------------------------------------------------------------------------\u001B[0m", - "\u001B[0;31mTypeError\u001B[0m Traceback (most recent call last)", - "Cell \u001B[0;32mIn[8], line 1\u001B[0m\n\u001B[0;32m----> 1\u001B[0m ctable \u001B[38;5;241m=\u001B[39m \u001B[43mCalibrationTable\u001B[49m\u001B[43m(\u001B[49m\u001B[43m)\u001B[49m\n", - "File \u001B[0;32m~/.conda/envs/calviper/lib/python3.10/site-packages/toolviper/utils/parameter.py:44\u001B[0m, in \u001B[0;36mvalidate.
<xarray.Dataset> Size: 26kB\n", + "Dimensions: (time: 20, antenna: 10, frequency: 8, polarization: 2,\n", + " gain: 1, scan_id: 20)\n", + "Coordinates:\n", + " * time (time) float64 160B 1.307e+09 1.307e+09 ... 1.307e+09\n", + " * antenna (antenna) <U9 360B 'DV02_A015' 'DV06_T704' ... 'PM03_J504'\n", + " * frequency (frequency) float64 64B 3.439e+11 3.439e+11 ... 3.44e+11\n", + " * polarization (polarization) <U1 8B 'X' 'Y'\n", + " * scan_id (scan_id) int64 160B 9 9 9 9 9 17 17 ... 21 21 25 25 25 25 25\n", + " * gain (gain) float64 8B 3.567e-319\n", + "Data variables:\n", + " GAIN (time, antenna, frequency, polarization, gain) float64 26kB ...\n", + "Attributes:\n", + " calibration_type: gain" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "execution_count": 9 }, { "cell_type": "code", diff --git a/src/calviper/factory/base.py b/src/calviper/factory/base.py new file mode 100644 index 0000000..49c73be --- /dev/null +++ b/src/calviper/factory/base.py @@ -0,0 +1,158 @@ +import xarray as xr + +import numpy as np + +from toolviper.utils import logger +from dataclasses import dataclass +from abc import ABC, abstractmethod +from typing import TypeVar, Type, Union + +T = TypeVar('T', bound='Parent') + + +@dataclass(init=False) +class PolarizationBasis: + name: str + type: str + length: int + + +class JonesMatrix(ABC): + + def __init__(self): + # private parent variables + #self._parameters: Union[np.array, None] = None + self._matrix: Union[np.array, None] = None + + # public parent variable + self.type: Union[str, None] = None + + #self.dtype: Union[type, None] = None + self.n_times: Union[int, None] = None + self.n_antennas: Union[int, None] = None + self.n_channels: Union[int, None] = None + self.n_polarizations: Union[int, None] = None + + #self.n_channel_matrices: Union[int, None] = None + self.n_parameters: Union[int, None] = None + self.caltable_name: Union[str, None] = None + + #self.channel_dependent_parameters: bool = False + + self.polarization_basis: PolarizationBasis = PolarizationBasis() + self.name: str = "BaseJonesMatrix" + + # Inherited member properties + @property + def shape(self) -> tuple: + raise NotImplementedError + + @shape.setter + def shape(self, shape: tuple): + return NotImplementedError + + @property + @abstractmethod + def parameters(self) -> np.ndarray: + raise NotImplementedError + + @parameters.setter + @abstractmethod + def parameters(self, array: np.array) -> None: + raise NotImplementedError + + @property + @abstractmethod + def matrix(self) -> np.ndarray: + return self._matrix + + @matrix.setter + @abstractmethod + def matrix(self, array: np.array) -> np.ndarray: + raise NotImplementedError + + @abstractmethod + def calculate(self) -> None: + raise NotImplementedError + + # Inherited method properties + @classmethod + def from_parameters(cls: Type[T], parameters: dict) -> T: + import inspect + + obj = cls() + updated_params = {} + + # This is a bit of a complicated way to do this BUT it should allow for a generic form of + # from_parameters() for all child classes. I think ... + for key, value in parameters.items(): + if key in inspect.getmembers(cls.__bases__[0], predicate=inspect.isfunction): + updated_params[f"_{key}"] = value + + elif key in inspect.getmembers(cls, predicate=inspect.isfunction): + updated_params[f"_{key}"] = value + + else: + if key in cls().__dict__.keys(): + updated_params[key] = value + + elif key in cls.__bases__[0]().__dict__.keys(): + updated_params[key] = value + + else: + pass + + vars(obj).update(updated_params) + + return obj + + @classmethod + def from_visibility(cls: Type[T], data: xr.Dataset, time_dependence: bool = False) -> T: + return cls + + # Commenting all these out for the moment. They were written in line with George's original code + # but that was based on a different workflow than I am foreseeing at the moment. These will be + # added back, if needed, as work progresses. + ''' + def initialize_parameters(self, dtype: np.dtype, shape: tuple = None): + # Set data type + self.type = dtype + + # Update shape is needed + if shape is not None: + self.shape = shape + + assert self.shape is not None, logger.error("Matrix shape is not set.") + + # Initialize the parameters to default + self.parameters = np.ones(self.shape, dtype=dtype) + + # Reset Jones + self.matrix = np.empty([]) + + def initialize_jones(self, shape: tuple = None): + if shape is not None: + self.shape = shape + + self.matrix = np.identity(2, dtype=np.complex64) + self.matrix = np.tile(self.matrix, [self.n_times, self.n_antennas, self.n_channel_matrices, 1, 1]) + + def invert(self) -> Union[np.array, None]: + if np.any(np.abs(np.linalg.det(self.matrix)) == 0.): + logger.error(f"Jones matrix is singular: {np.linalg.det(self.matrix)}") + return None + + return np.linalg.inv(self.matrix) + + def accumulate(self, other: Type[T]) -> T: + # I think this could just be an overload of __mul__() + return np.matmul(other.matrix, self.matrix, out=self.matrix) + + def apply_left(self): + # Need to inspect use case + pass + + def apply_right(self): + # Need to inspect use case + pass +''' \ No newline at end of file diff --git a/src/calviper/factory/jones.py b/src/calviper/factory/jones.py new file mode 100644 index 0000000..ebc3cec --- /dev/null +++ b/src/calviper/factory/jones.py @@ -0,0 +1,103 @@ +import numpy as np +import xarray as xr + +import calviper.math.tools as tools +import toolviper.utils.logger as logger + +from calviper.base import JonesMatrix +from typing import TypeVar, Type, Union + +T = TypeVar('T', bound='Parent') + + +class GainJones(JonesMatrix): + def __init__(self): + super(GainJones, self).__init__() + + # Public parent variable + self.n_times = None + self.type: Union[str, None] = "G" + + #self.dtype = np.complex64 + self.n_polarizations: Union[int, None] = 4 + self.n_parameters: Union[int, None] = None + self.n_baselines: Union[int, None] = None + self.n_channels: Union[int, None] = None + self.channel_dependent_parameters: bool = False + + # Private variables + self._parameters = None + self._matrix = None + self._antenna_map = None + + self.name: str = "GainJonesMatrix" + + # This is just an example of how this would be done. There should certainly be checks and customization + # but for now just set the values simply as the original code doesn't do anything more complicated for now. + @property + def parameters(self) -> np.ndarray: + return self._parameters + + @parameters.setter + def parameters(self, array: np.ndarray) -> None: + self._parameters = array + + @property + def matrix(self) -> np.ndarray: + return self._matrix + + @matrix.setter + def matrix(self, array: np.ndarray) -> None: + #(self.n_times, self.n_baselines, self.n_channels, _, _) = array.shape + + # There should be a check on the shape here. I don't think we want to allow, for instance, + # an axis to be averaged while also having the dimensions stored in the object not change. + self._matrix = array + + def calculate(self) -> None: + #self.initialize_jones() + + self.matrix = np.identity(2, dtype=np.complex64) + self.matrix = np.tile(self.matrix, [self.n_times, self.n_antennas, self.n_channel_matrices, 1, 1]) + + @classmethod + def from_visibility(cls: Type[T], dataset: xr.Dataset, time_dependence: bool = False) -> T: + """ + Build Jones matrix from visibility data. + :param dataset: + :param time_dependence: + :return: + """ + + shape = dataset.VISIBILITY.shape + + # This will encode the antenna values into an integer list. + index, antennas = tools.encode(dataset.baseline_antenna1_name.to_numpy()) + + instance = cls() + + # There should be a gain value for each independent antenna. Here we choose antenna_1 names but either + # would work fine. + instance.n_antennas = np.unique(dataset.baseline_antenna1_name).shape[0] + + # With no polarization and one channel, n_parameters = n_antennas + # instance.n_parameters = n_parameters + instance.n_parameters = instance.n_antennas * instance.n_polarizations + + polarization_axis_ = int(instance.n_polarizations // 2) + + identity = np.identity(polarization_axis_, dtype=np.complex64) + + instance._antenna_map = {i: str(antenna) for i, antenna in enumerate(antennas)} + + instance.n_times, instance.n_baselines, instance.n_channels, instance.n_polarizations = shape + + # Initialize default parameters + instance.parameters = np.empty((instance.n_times, instance.n_channels, instance.n_parameters), + dtype=np.complex64) + + # Build on the above idea ... wrong as they may be. Simplicity first. + # instance.matrix = np.tile(identity, reps=[*shape, 1, 1]) + instance.matrix = np.tile(identity, reps=[instance.n_times, instance.n_baselines, instance.n_channels, 1, 1]) + + return instance diff --git a/src/calviper/factory/table.py b/src/calviper/factory/table.py new file mode 100644 index 0000000..6b6595b --- /dev/null +++ b/src/calviper/factory/table.py @@ -0,0 +1,101 @@ +import numpy as np +import xarray as xr + +import toolviper.utils.logger as logger +import toolviper.utils.parameter + +from abc import ABC +from abc import abstractmethod + +from typing import Union + + +class BaseJonesMatrix(ABC): + + # Base calibration table abstract class + @abstractmethod + def generate(self, coords: dict) -> Union[xr.Dataset, None]: + pass + + +class JonesFactory(ABC): + # Base factory class for table factory + @abstractmethod + def create_table(self, factory: Union[None, str]): + pass + + +class GainTable(BaseJonesMatrix): + + # This is intended to be an implementation of a gain table simulator. It is + # currently very rough and filled with random numbers. Generally based on the + # original cal.py + def generate(self, coords: dict) -> xr.Dataset: + shape = tuple(value.shape[0] for value in coords.values()) + + dims = {} + for key, value in coords.items(): + dims[key] = value.shape[0] + + parameter = np.random.uniform(-np.pi, np.pi, shape) + amplitude = np.random.normal(1.0, 0.1, shape) + parameter = np.vectorize(complex)( + np.cos(parameter), + np.sin(parameter) + ) + + xds = xr.Dataset() + + xds["PARAMETER"] = xr.DataArray(amplitude * parameter, dims=dims) + xds = xds.assign_coords(coords) + + return xds + + @staticmethod + def empty_like(dataset: xr.Dataset) -> xr.Dataset: + antenna = dataset.antenna_xds.antenna_name.values + polarizations = np.unique([p for value in dataset.polarization.values for p in list(value)]) + + dims = dict( + time=dataset.sizes["time"], + antenna=antenna.shape[0], + frequency=dataset.sizes["frequency"], + polarization=polarizations.shape[0], + gain=1 + ) + + coords = dict( + time=(["time"], dataset.time.values), + antenna=(["antenna"], antenna), + frequency=(["frequency"], dataset.frequency.values), + polarization=(["polarization"], polarizations), + scan_id=(["scan_id"], dataset.scan_number.values), + gain=(["gain"], np.empty(1)) + ) + + gain = np.empty(list(dims.values())) + + xds = xr.Dataset() + + xds["PARAMETER"] = xr.DataArray(gain, dims=dims) + xds.attrs["calibration_type"] = "gain" + xds = xds.assign_coords(coords) + + return xds + + +class CalibrationTable(JonesFactory): + + def __init__(self): + self.factory_list = { + "gain": GainTable, + } + + @toolviper.utils.parameter.validate() + def create_table(self, factory: str) -> Union[BaseJonesMatrix, None]: + try: + return self.factory_list[factory]() + + except KeyError: + logger.error(f"Factory method, {factory} not implemented.") + return None