From d2403a84a66da0d6b0c5acddbe8b0eff217c4178 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 5 Aug 2021 15:38:27 +0200 Subject: [PATCH 01/89] Data import: potentially ignore flags --- misc/ms2npz.py | 2 ++ resolve/data/ms_import.py | 24 +++++++++++++++++------- 2 files changed, 19 insertions(+), 7 deletions(-) diff --git a/misc/ms2npz.py b/misc/ms2npz.py index a9e6ee77..598344a4 100644 --- a/misc/ms2npz.py +++ b/misc/ms2npz.py @@ -18,6 +18,7 @@ if __name__ == "__main__": parser.add_argument("--ch-begin", type=int) parser.add_argument("--ch-end", type=int) parser.add_argument("--ch-jump", type=int) + parser.add_argument("--ignore-flags", action="store_true") parser.add_argument("ms") parser.add_argument( "polarization_mode", @@ -35,6 +36,7 @@ if __name__ == "__main__": args.spectral_window, args.polarization_mode, slice(args.ch_begin, args.ch_end, args.ch_jump), + args.ignore_flags ) for ifield, oo in enumerate(obs): fname = join(args.output_folder, f"{name}field{ifield}.npz") diff --git a/resolve/data/ms_import.py b/resolve/data/ms_import.py index 34105deb..bef17c8b 100644 --- a/resolve/data/ms_import.py +++ b/resolve/data/ms_import.py @@ -24,6 +24,7 @@ def ms2observations( spectral_window, polarizations="all", channel_slice=slice(None), + ignore_flags=False, ): """Read and convert a given measurement set into an array of :class:`Observation` @@ -50,6 +51,9 @@ def ms2observations( channel_slice : slice Slice of selected channels. Default: select all channels FIXME Select channels by indices + ignore_flags : bool + If True, the whole measurement set is imported irrespective of the + flags. Default is false. Returns ------- @@ -131,6 +135,7 @@ def ms2observations( pol_summation, with_calib_info, channel_slice, + ignore_flags, ) vis[wgt == 0] = 0.0 vis = _ms2resolve_transpose(vis) @@ -167,6 +172,7 @@ def read_ms_i( pol_summation, with_calib_info, channel_slice, + ignore_flags, ): freq = np.array(freq) my_asserteq(freq.ndim, 1) @@ -197,9 +203,7 @@ def read_ms_i( while start < nrow: print("First pass:", f"{(start/nrow*100):.1f}%", end="\r") stop = min(nrow, start + step) - tflags = t.getcol("FLAG", startrow=start, nrow=stop - start)[ - ..., pol_indices - ] + tflags = _conditional_flags(t, start, stop, pol_indices, ignore_flags) twgt = t.getcol(weightcol, startrow=start, nrow=stop - start)[ ..., pol_indices ] @@ -273,16 +277,15 @@ def read_ms_i( tvis = tvis[active_rows[start:stop]] tvis = tvis[:, active_channels] - tflags = t.getcol("FLAG", startrow=start, nrow=stop - start)[ - ..., pol_indices - ] + tflags = _conditional_flags(t, start, stop, pol_indices, ignore_flags) if not allrows: tflags = tflags[active_rows[start:stop]] tflags = tflags[:, active_channels] assert twgt.ndim == tflags.ndim == 3 assert tflags.dtype == np.bool - twgt = twgt * (~tflags) + if not ignore_flags: + twgt = twgt * (~tflags) if pol_summation: assert twgt.shape[2] == 2 # Noise-weighted average @@ -322,3 +325,10 @@ def ms_n_spectral_windows(ms): with ms_table(join(ms, "SPECTRAL_WINDOW")) as t: freq = t.getcol("CHAN_FREQ") return freq.shape[0] + + +def _conditional_flags(table, start, stop, pol_indices, ignore): + tflags = table.getcol("FLAG", startrow=start, nrow=stop - start)[..., pol_indices] + if ignore: + tflags = np.zeros_like(tflags) + return tflags -- GitLab From e6505e23fdb2439fe0d6a36507fb478961e05e85 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 6 Aug 2021 17:12:08 +0200 Subject: [PATCH 02/89] Support export of autocorrelations only --- misc/ms2npz.py | 10 +++++++++- resolve/data/observation.py | 4 ++++ 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/misc/ms2npz.py b/misc/ms2npz.py index 598344a4..811cc287 100644 --- a/misc/ms2npz.py +++ b/misc/ms2npz.py @@ -19,6 +19,7 @@ if __name__ == "__main__": parser.add_argument("--ch-end", type=int) parser.add_argument("--ch-jump", type=int) parser.add_argument("--ignore-flags", action="store_true") + parser.add_argument("--autocorrelations-only", action="store_true") parser.add_argument("ms") parser.add_argument( "polarization_mode", @@ -39,6 +40,13 @@ if __name__ == "__main__": args.ignore_flags ) for ifield, oo in enumerate(obs): - fname = join(args.output_folder, f"{name}field{ifield}.npz") + if args.autocorrelations_only: + oo = oo.restrict_to_autocorrelations() + auto = "autocorrelationsonly" + else: + auto = "" + fname = join(args.output_folder, f"{name}field{ifield}{args.data_column}{auto}.npz") print(f"Save {fname}") oo.save(fname, args.compress) + if oo.vis.size == 0: + print(f"WARNING: {fname} is empty") diff --git a/resolve/data/observation.py b/resolve/data/observation.py index 2dbe6c15..ab2e16fd 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -433,6 +433,10 @@ class Observation(BaseObservation): pol = self._polarization.restrict_to_stokes_i() return Observation(self._antpos, vis, wgt, pol, self._freq, self._direction) + def restrict_to_autocorrelations(self): + slc = self._antpos.ant1 == self._antpos.ant2 + return self[slc] + def move_time(self, t0): antpos = self._antpos.move_time(t0) return Observation( -- GitLab From 102bf40b1ba859afc7b8202617f054d63864bfdd Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Tue, 10 Aug 2021 12:31:12 +0200 Subject: [PATCH 03/89] Add documentation --- resolve/multi_frequency/operators.py | 33 +++++++++++++++++----------- 1 file changed, 20 insertions(+), 13 deletions(-) diff --git a/resolve/multi_frequency/operators.py b/resolve/multi_frequency/operators.py index 583e6707..067d00ba 100644 --- a/resolve/multi_frequency/operators.py +++ b/resolve/multi_frequency/operators.py @@ -1,29 +1,36 @@ # SPDX-License-Identifier: GPL-3.0-or-later -# Copyright(C) 2013-2020 Max-Planck-Society +# Copyright(C) 2013-2021 Max-Planck-Society # Authors: Philipp Frank, Philipp Arras, Philipp Haim import numpy as np import nifty8 as ift -from ..util import my_asserteq +from ..util import my_assert, my_assert_isinstance, my_asserteq +from .irg_space import IRGSpace class WienerIntegrations(ift.LinearOperator): - def __init__(self, freqdomain, imagedomain): - # FIXME Write interface checks - self._target = ift.makeDomain((freqdomain, imagedomain)) - dom = list(self._target) - dom = ift.UnstructuredDomain((2, freqdomain.size - 1)), imagedomain + """Operator that performs the integrations necessary for an integrated + Wiener process. + + Parameters + ---------- + time_domain : IRGSpace + Domain that contains the temporal information of the process. + + remaining_domain : DomainTuple or Domain + All integrations are handled independently for this domain. + """ + def __init__(self, time_domain, remaining_domain): + my_assert_isinstance(time_domain, IRGSpace) + self._target = ift.makeDomain((time_domain, remaining_domain)) + dom = ift.UnstructuredDomain((2, time_domain.size - 1)), remaining_domain self._domain = ift.makeDomain(dom) - self._volumes = freqdomain.dvol[:, None, None] + self._volumes = time_domain.dvol[:, None, None] self._capability = self.TIMES | self.ADJOINT_TIMES def apply(self, x, mode): - # FIXME If it turns out that this operator is a performance - # bottleneck we can try implement it in parallel in C++. But - # it may be hard to achieve good scaling because I think it - # becomes memory bound quickly. self._check_input(x, mode) first, second = (0,), (1,) from_second = (slice(1, None),) @@ -50,7 +57,7 @@ class WienerIntegrations(ift.LinearOperator): def IntWProcessInitialConditions(a0, b0, wfop): for op in [a0, b0, wfop]: - ift.is_operator(op) + my_assert(ift.is_operator(op)) my_asserteq(a0.target, b0.target, ift.makeDomain(wfop.target[1])) bc = _FancyBroadcast(wfop.target) factors = ift.full(wfop.target[0], 0) -- GitLab From 2bf1fda1b30d02b17195ee3fcc4a41c19ab831e5 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Tue, 10 Aug 2021 15:42:55 +0200 Subject: [PATCH 04/89] MfWeightingInterpolation takes Field --- resolve/multi_frequency/operators.py | 2 +- test/test_general.py | 8 +++++--- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/resolve/multi_frequency/operators.py b/resolve/multi_frequency/operators.py index 067d00ba..11aae3db 100644 --- a/resolve/multi_frequency/operators.py +++ b/resolve/multi_frequency/operators.py @@ -97,7 +97,7 @@ class MfWeightingInterpolation(ift.LinearOperator): # FIXME Try to unify all those operators which loop over freq dimension self._ops = [] for ii in range(nfreq): - op = ift.LinearInterpolator(domain[1], eff_uvw[:, :, ii]) + op = ift.LinearInterpolator(domain[1], eff_uvw.val[:, :, ii]) self._ops.append(op) my_asserteq(self.target.shape[0], 1) diff --git a/test/test_general.py b/test/test_general.py index ec799a03..e5e2f997 100644 --- a/test/test_general.py +++ b/test/test_general.py @@ -309,9 +309,11 @@ def test_randomonmaster(): def test_mfweighting(): nrow = 100 nchan = 4 - effuv = ift.random.current_rng().random((nrow, nchan)) - dom = ift.UnstructuredDomain(nchan), ift.RGSpace(npix, 2 * np.max(effuv) / npix) - op = rve.MfWeightingInterpolation(effuv[None], dom) + npol = 1 + effuv = ift.random.current_rng().random((npol, nrow, nchan)) + effuv = ift.makeField([ift.UnstructuredDomain(ii) for ii in effuv.shape], effuv) + dom = ift.UnstructuredDomain(nchan), ift.RGSpace(npix, 2 * np.max(effuv.val) / npix) + op = rve.MfWeightingInterpolation(effuv, dom) ift.extra.check_linear_operator(op) -- GitLab From 82503ac179bed53f094122d0ba85d6dba9cd5e8e Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 12 Aug 2021 12:07:44 +0200 Subject: [PATCH 05/89] Cosmetics --- misc/ms2npz.py | 1 + resolve/data/ms_import.py | 29 +++++++---------------------- 2 files changed, 8 insertions(+), 22 deletions(-) diff --git a/misc/ms2npz.py b/misc/ms2npz.py index 811cc287..cbe817c8 100644 --- a/misc/ms2npz.py +++ b/misc/ms2npz.py @@ -1,3 +1,4 @@ +#!/usr/bin/python3 # SPDX-License-Identifier: GPL-3.0-or-later # Copyright(C) 2019-2020 Max-Planck-Society # Author: Philipp Arras diff --git a/resolve/data/ms_import.py b/resolve/data/ms_import.py index bef17c8b..3bf5beca 100644 --- a/resolve/data/ms_import.py +++ b/resolve/data/ms_import.py @@ -1,15 +1,16 @@ # SPDX-License-Identifier: GPL-3.0-or-later # Copyright(C) 2019-2021 Max-Planck-Society +import os from os.path import isdir, join, splitext import numpy as np +from ..util import my_assert, my_assert_isinstance, my_asserteq from .antenna_positions import AntennaPositions from .direction import Direction from .observation import Observation from .polarization import Polarization -from ..util import my_assert, my_asserteq, my_assert_isinstance def ms_table(path): @@ -17,15 +18,9 @@ def ms_table(path): return table(path, readonly=True, ack=False) -def ms2observations( - ms, - data_column, - with_calib_info, - spectral_window, - polarizations="all", - channel_slice=slice(None), - ignore_flags=False, -): +def ms2observations(ms, data_column, with_calib_info, spectral_window, + polarizations="all", channel_slice=slice(None), + ignore_flags=False): """Read and convert a given measurement set into an array of :class:`Observation` If WEIGHT_SPECTRUM is available this column is used for weighting. @@ -162,18 +157,8 @@ def _determine_weighting(t): return fullwgt, weightcol -def read_ms_i( - name, - data_column, - freq, - field, - spectral_window, - pol_indices, - pol_summation, - with_calib_info, - channel_slice, - ignore_flags, -): +def read_ms_i(name, data_column, freq, field, spectral_window, pol_indices, + pol_summation, with_calib_info, channel_slice, ignore_flags): freq = np.array(freq) my_asserteq(freq.ndim, 1) my_assert(len(freq) > 0) -- GitLab From 27dce8a0a1faa885acdbe8642cccfe970becbb86 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 12 Aug 2021 12:07:55 +0200 Subject: [PATCH 06/89] Support "." as ms location --- resolve/data/ms_import.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/resolve/data/ms_import.py b/resolve/data/ms_import.py index 3bf5beca..6cb27daa 100644 --- a/resolve/data/ms_import.py +++ b/resolve/data/ms_import.py @@ -63,7 +63,11 @@ def ms2observations(ms, data_column, with_calib_info, spectral_window, """ if ms[-1] == "/": ms = ms[:-1] - if not isdir(ms) or splitext(ms)[1].lower() != ".ms": + if not isdir(ms): + raise RuntimeError + if ms == ".": + ms = os.getcwd() + if splitext(ms)[1].lower() != ".ms": raise RuntimeError if isinstance(polarizations, str): polarizations = [polarizations] -- GitLab From 6d789a3e346f91061a3b96625476ac189362e1fe Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 12 Aug 2021 12:08:10 +0200 Subject: [PATCH 07/89] Docs --- misc/ms2npz.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/misc/ms2npz.py b/misc/ms2npz.py index cbe817c8..97ee29b5 100644 --- a/misc/ms2npz.py +++ b/misc/ms2npz.py @@ -20,7 +20,9 @@ if __name__ == "__main__": parser.add_argument("--ch-end", type=int) parser.add_argument("--ch-jump", type=int) parser.add_argument("--ignore-flags", action="store_true") - parser.add_argument("--autocorrelations-only", action="store_true") + parser.add_argument("--autocorrelations-only", action="store_true", + help=("If this flag is set, all autocorrelations are imported " + "irrespective of whether they are flagged or not.")) parser.add_argument("ms") parser.add_argument( "polarization_mode", -- GitLab From 831bbffde5f165c0e10a5c0bfc225a0d38c3b7c9 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 12 Aug 2021 13:14:57 +0200 Subject: [PATCH 08/89] Add support for auxiliary MS tables --- misc/ms2npz.py | 0 resolve/data/auxiliary_table.py | 45 +++++++++++++++++++++++++++++++ resolve/data/ms_import.py | 11 +++++++- resolve/data/observation.py | 47 +++++++++++++++++++++++++++------ 4 files changed, 94 insertions(+), 9 deletions(-) mode change 100644 => 100755 misc/ms2npz.py create mode 100644 resolve/data/auxiliary_table.py diff --git a/misc/ms2npz.py b/misc/ms2npz.py old mode 100644 new mode 100755 diff --git a/resolve/data/auxiliary_table.py b/resolve/data/auxiliary_table.py new file mode 100644 index 00000000..eb06aea9 --- /dev/null +++ b/resolve/data/auxiliary_table.py @@ -0,0 +1,45 @@ +# SPDX-License-Identifier: GPL-3.0-or-later +# Copyright(C) 2021 Max-Planck-Society +# Author: Philipp Arras + +import numpy as np + +from ..util import compare_attributes, my_assert, my_assert_isinstance, my_asserteq + + +class AuxiliaryTable: + def __init__(self, inp_dict): + my_assert_isinstance(inp_dict, dict) + nrow = None + for kk, lst in inp_dict.items(): + my_assert_isinstance(kk, str) + my_assert_isinstance(lst, (list, np.ndarray)) + if nrow is None: + nrow = len(lst) + my_asserteq(nrow, len(lst)) + my_asserteq(type(elem) for elem in lst) + self._dct = inp_dict + self._eq_attributes = [self._dct] + + def __eq__(self, other): + if not isinstance(other, type(self)): + return False + return compare_attributes(self, other, self._eq_attributes) + + def __str__(self): + s = ["AuxiliaryTable:"] + for kk, lst in self._dct.items(): + if isinstance(lst, np.ndarray): + s.append(f" {kk:<20} {str(lst.shape):>10} {str(lst.dtype):>15}") + else: + s.append(f" {kk:<20} {len(lst):>10} {str(type(lst[0])):>15}") + return "\n".join(s) + + def to_list(self): + return [list(self._dct.keys())] + [ss for ss in self._dct.values()] + + @staticmethod + def from_list(lst): + keys = lst[0] + dct = {kk: lst[ii+1] for ii, kk in enumerate(keys)} + return AuxiliaryTable(dct) diff --git a/resolve/data/ms_import.py b/resolve/data/ms_import.py index 6cb27daa..b0a48ed7 100644 --- a/resolve/data/ms_import.py +++ b/resolve/data/ms_import.py @@ -8,6 +8,7 @@ import numpy as np from ..util import my_assert, my_assert_isinstance, my_asserteq from .antenna_positions import AntennaPositions +from .auxiliary_table import AuxiliaryTable from .direction import Direction from .observation import Observation from .polarization import Polarization @@ -120,6 +121,12 @@ def ms2observations(ms, data_column, with_calib_info, spectral_window, dirs.append(Direction(pc[0], equinox)) dirs = tuple(dirs) + auxtables = {} + with ms_table(join(ms, "ANTENNA")) as t: + keys = ["NAME", "STATION", "TYPE", "MOUNT", "POSITION", "OFFSET", + "DISH_DIAMETER"] + auxtables["ANTENNA"] = AuxiliaryTable({kk: t.getcol(kk) for kk in keys}) + # FIXME Determine which observation is calibration observation # FIXME Import name of source observations = [] @@ -140,7 +147,9 @@ def ms2observations(ms, data_column, with_calib_info, spectral_window, vis = _ms2resolve_transpose(vis) wgt = _ms2resolve_transpose(wgt) antpos = AntennaPositions(uvw, ant1, ant2, time) - observations.append(Observation(antpos, vis, wgt, polobj, freq_out, direction)) + obs = Observation(antpos, vis, wgt, polobj, freq_out, direction, + auxiliary_tables=auxtables) + observations.append(obs) return observations diff --git a/resolve/data/observation.py b/resolve/data/observation.py index ab2e16fd..e0455ae6 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -2,16 +2,17 @@ # Copyright(C) 2019-2021 Max-Planck-Society # Author: Philipp Arras -import numpy as np - import nifty8 as ift +import numpy as np -from .antenna_positions import AntennaPositions from ..constants import SPEEDOFLIGHT -from .direction import Direction, Directions from ..mpi import onlymaster +from ..util import (compare_attributes, my_assert, my_assert_isinstance, + my_asserteq) +from .antenna_positions import AntennaPositions +from .auxiliary_table import AuxiliaryTable +from .direction import Direction, Directions from .polarization import Polarization -from ..util import compare_attributes, my_assert, my_assert_isinstance, my_asserteq class BaseObservation: @@ -157,7 +158,7 @@ class BaseObservation: return NotImplementedError def __eq__(self, other): - if not isinstance(other, self.__class__): + if not isinstance(other, type(self)): return False if ( self._vis.dtype != other._vis.dtype @@ -268,13 +269,15 @@ class Observation(BaseObservation): Contains the measured frequencies. Shape (n_channels) direction : Direction Direction information of the data set. + auxiliary_tables : dict + Dictionary of auxiliary tables. Default: None. Note ---- vis and weight must have the same shape. """ - def __init__(self, antenna_positions, vis, weight, polarization, freq, direction): + def __init__(self, antenna_positions, vis, weight, polarization, freq, direction, auxiliary_tables=None): nrows = len(antenna_positions) my_assert_isinstance(direction, Direction) my_assert_isinstance(polarization, Polarization) @@ -304,8 +307,15 @@ class Observation(BaseObservation): self._freq = freq self._direction = direction - self._eq_attributes = "_direction", "_polarization", "_freq", "_antpos", "_vis", "_weight" + if auxiliary_tables is not None: + my_assert_isinstance(auxiliary_tables, dict) + for kk, vv in auxiliary_tables.items(): + my_assert_isinstance(vv, AuxiliaryTable) + my_assert_isinstance(kk, str) + self._auxiliary_tables = auxiliary_tables + self._eq_attributes = ("_direction", "_polarization", "_freq", + "_antpos", "_vis", "_weight", "_auxiliary_tables") @onlymaster def save(self, file_name, compress): @@ -320,6 +330,10 @@ class Observation(BaseObservation): if vv is None: vv = np.array([]) dct[f"antpos{ii}"] = vv + if self._auxiliary_tables is not None: + for kk, auxtable in self._auxiliary_tables.items(): + for ii, elem in enumerate(auxtable.to_list()): + dct[f"auxtable_{kk}_{ii}"] = elem f = np.savez_compressed if compress else np.savez f(file_name, **dct) @@ -332,6 +346,22 @@ class Observation(BaseObservation): if val.size == 0: val = None antpos.append(val) + + # Load auxtables + keys = set(kk.split("_")[1] for kk in dct.keys() if kk[:8] == "auxtable") + if len(keys) == 0: + auxtables = None + else: + auxtables = {} + for kk in keys: + ii = 0 + inp = [] + while f"auxtable_{kk}_{ii}" in dct.keys(): + inp.append(dct[f"auxtable_{kk}_{ii}"]) + ii += 1 + auxtables[kk] = AuxiliaryTable.from_list(inp) + # /Load auxtables + pol = Polarization.from_list(dct["polarization"]) direction = Direction.from_list(dct["direction"]) slc = slice(None) if lo_hi_index is None else slice(*lo_hi_index) @@ -343,6 +373,7 @@ class Observation(BaseObservation): pol, dct["freq"][slc], direction, + auxtables ) def flags_to_nan(self): -- GitLab From a46de75524b4bc8592ccb36cb75dc793d56f246f Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 12 Aug 2021 15:30:40 +0200 Subject: [PATCH 09/89] Improve comparison --- resolve/data/auxiliary_table.py | 2 +- resolve/util.py | 54 ++++++++++++++++++++++++--------- 2 files changed, 41 insertions(+), 15 deletions(-) diff --git a/resolve/data/auxiliary_table.py b/resolve/data/auxiliary_table.py index eb06aea9..58365b8c 100644 --- a/resolve/data/auxiliary_table.py +++ b/resolve/data/auxiliary_table.py @@ -19,7 +19,7 @@ class AuxiliaryTable: my_asserteq(nrow, len(lst)) my_asserteq(type(elem) for elem in lst) self._dct = inp_dict - self._eq_attributes = [self._dct] + self._eq_attributes = "_dct", def __eq__(self, other): if not isinstance(other, type(self)): diff --git a/resolve/util.py b/resolve/util.py index 4c344455..b8eda2b3 100644 --- a/resolve/util.py +++ b/resolve/util.py @@ -28,20 +28,46 @@ def my_assert_isinstance(*args): def compare_attributes(obj0, obj1, attribute_list): - for a in attribute_list: - compare = getattr(obj0, a) != getattr(obj1, a) - if type(compare) is bool: - if compare: - return False - continue - if isinstance(compare, ift.MultiField): - for vv in compare.val.values(): - if vv.any(): - return False - continue - if compare.any(): - return False - return True + return all(_fancy_equal(getattr(obj0, a), getattr(obj1, a)) + for a in attribute_list) + + +def _fancy_equal(o1, o2): + if not _types_equal(o1, o2): + return False + + # Turn MultiField into dict + if isinstance(o1, ift.MultiField): + o1, o2 = o1.val, o2.val + + # Compare dicts + if isinstance(o1, dict): + return _deep_equal(o1, o2) + + # Compare simple objects and np.ndarrays + return _compare_simple_or_array(o1, o2) + + +def _deep_equal(a, b): + if not isinstance(a, dict) or not isinstance(b, dict): + raise TypeError + + if a.keys() != b.keys(): + return False + + return all(_compare_simple_or_array(a[kk], b[kk]) for kk in a.keys()) + + +def _compare_simple_or_array(a, b): + equal = a == b + try: + return bool(equal) + except ValueError: + return equal.all() + + +def _types_equal(a, b): + return type(a) == type(b) class Reshaper(ift.LinearOperator): -- GitLab From 64dfd86c4dba93cccc420aa5e698b5301f6b1404 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 12 Aug 2021 15:35:56 +0200 Subject: [PATCH 10/89] Expose auxiliary tables --- resolve/data/observation.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/resolve/data/observation.py b/resolve/data/observation.py index e0455ae6..a3a561e2 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -487,6 +487,9 @@ class Observation(BaseObservation): def antenna_positions(self): return self._antpos + def auxiliary_table(name): + return self._auxiliary_tables[name] + def effective_uvw(self): out = np.einsum("ij,k->jik", self.uvw, self._freq / SPEEDOFLIGHT) my_asserteq(out.shape, (3, self.nrow, self.nfreq)) -- GitLab From 55bc2826d8d7e33b8fa7de4849e8dc069e7c8e85 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 12 Aug 2021 18:28:36 +0200 Subject: [PATCH 11/89] Add more aux tables --- resolve/data/ms_import.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/resolve/data/ms_import.py b/resolve/data/ms_import.py index b0a48ed7..292f7b69 100644 --- a/resolve/data/ms_import.py +++ b/resolve/data/ms_import.py @@ -123,9 +123,17 @@ def ms2observations(ms, data_column, with_calib_info, spectral_window, auxtables = {} with ms_table(join(ms, "ANTENNA")) as t: - keys = ["NAME", "STATION", "TYPE", "MOUNT", "POSITION", "OFFSET", - "DISH_DIAMETER"] + keys = ["NAME", "STATION", "TYPE", "MOUNT", "POSITION", "OFFSET", "DISH_DIAMETER"] auxtables["ANTENNA"] = AuxiliaryTable({kk: t.getcol(kk) for kk in keys}) + with ms_table(join(ms, "FIELD")) as t: + keys = ["NAME", "CODE", "TIME", "NUM_POLY", "DELAY_DIR", "PHASE_DIR", "REFERENCE_DIR", + "SOURCE_ID"] + auxtables["FIELD"] = AuxiliaryTable({kk: t.getcol(kk) for kk in keys}) + with ms_table(join(ms, "SPECTRAL_WINDOW")) as t: + keys = ["NAME", "REF_FREQUENCY", "CHAN_FREQ", "CHAN_WIDTH", "MEAS_FREQ_REF", "EFFECTIVE_BW", + "RESOLUTION", "TOTAL_BANDWIDTH", "NET_SIDEBAND", "IF_CONV_CHAIN", "FREQ_GROUP", + "FREQ_GROUP_NAME"] + auxtables["SPECTRAL_WINDOW"] = AuxiliaryTable({kk: t.getcol(kk) for kk in keys}) # FIXME Determine which observation is calibration observation # FIXME Import name of source -- GitLab From 21fcc98bfebdd2c8e121cb75014e8a85927aa84d Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 13 Aug 2021 09:09:18 +0200 Subject: [PATCH 12/89] fixup! Add more aux tables --- resolve/data/observation.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/resolve/data/observation.py b/resolve/data/observation.py index a3a561e2..fa602c42 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -333,7 +333,7 @@ class Observation(BaseObservation): if self._auxiliary_tables is not None: for kk, auxtable in self._auxiliary_tables.items(): for ii, elem in enumerate(auxtable.to_list()): - dct[f"auxtable_{kk}_{ii}"] = elem + dct[f"auxtable_{kk}_{ii:>04}"] = elem f = np.savez_compressed if compress else np.savez f(file_name, **dct) @@ -348,7 +348,7 @@ class Observation(BaseObservation): antpos.append(val) # Load auxtables - keys = set(kk.split("_")[1] for kk in dct.keys() if kk[:8] == "auxtable") + keys = set(kk[9:-5] for kk in dct.keys() if kk[:9] == "auxtable_") if len(keys) == 0: auxtables = None else: @@ -356,8 +356,8 @@ class Observation(BaseObservation): for kk in keys: ii = 0 inp = [] - while f"auxtable_{kk}_{ii}" in dct.keys(): - inp.append(dct[f"auxtable_{kk}_{ii}"]) + while f"auxtable_{kk}_{ii:>04}" in dct.keys(): + inp.append(dct[f"auxtable_{kk}_{ii:>04}"]) ii += 1 auxtables[kk] = AuxiliaryTable.from_list(inp) # /Load auxtables -- GitLab From 51e6fbb5f8ef223877c359a9208f8d4dad846114 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 13 Aug 2021 09:12:29 +0200 Subject: [PATCH 13/89] Fixup --- resolve/data/observation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/resolve/data/observation.py b/resolve/data/observation.py index fa602c42..8e035487 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -487,7 +487,7 @@ class Observation(BaseObservation): def antenna_positions(self): return self._antpos - def auxiliary_table(name): + def auxiliary_table(self, name): return self._auxiliary_tables[name] def effective_uvw(self): -- GitLab From cd2bf383614517fcd39afb33c69b86b0d3378077 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 13 Aug 2021 10:32:05 +0200 Subject: [PATCH 14/89] More general ms_analysis --- misc/ms_summary.py | 12 ++++++++++++ resolve/data/ms_import.py | 4 ++-- 2 files changed, 14 insertions(+), 2 deletions(-) mode change 100644 => 100755 misc/ms_summary.py diff --git a/misc/ms_summary.py b/misc/ms_summary.py old mode 100644 new mode 100755 index c36bed83..20f781c0 --- a/misc/ms_summary.py +++ b/misc/ms_summary.py @@ -1,8 +1,10 @@ +#!/usr/bin/python3 # SPDX-License-Identifier: GPL-3.0-or-later # Copyright(C) 2021 Max-Planck-Society # Author: Philipp Arras import argparse +from os.path import join import resolve as rve @@ -16,3 +18,13 @@ if __name__ == "__main__": print() print(f"The data set has {nspec} spectral windows.") print() + + with rve.ms_table(join(args.ms, "FIELD")) as t: + name = t.getcol("NAME") + refdir = t.getcol("REFERENCE_DIR") + deldir = t.getcol("DELAY_DIR") + phdir = t.getcol("PHASE_DIR") + + print("NAME REFERNCE_DIR DELAY_DIR PHASE_DIR") + for nn, rd, dd, pd in zip(name, refdir, deldir, phdir): + print(nn, rd, dd, pd) diff --git a/resolve/data/ms_import.py b/resolve/data/ms_import.py index 292f7b69..aa7fe428 100644 --- a/resolve/data/ms_import.py +++ b/resolve/data/ms_import.py @@ -329,8 +329,8 @@ def read_ms_i(name, data_column, freq, field, spectral_window, pol_indices, def ms_n_spectral_windows(ms): with ms_table(join(ms, "SPECTRAL_WINDOW")) as t: - freq = t.getcol("CHAN_FREQ") - return freq.shape[0] + n_spectral_windows = t.nrows() + return n_spectral_windows def _conditional_flags(table, start, stop, pol_indices, ignore): -- GitLab From 341f972760a55907d343a6980fc6cafb462e897d Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 13 Aug 2021 11:32:19 +0200 Subject: [PATCH 15/89] Cosmetics --- misc/ms2npz.py | 5 +++-- resolve/data/ms_import.py | 17 +++++------------ 2 files changed, 8 insertions(+), 14 deletions(-) diff --git a/misc/ms2npz.py b/misc/ms2npz.py index 97ee29b5..7c4fa716 100755 --- a/misc/ms2npz.py +++ b/misc/ms2npz.py @@ -21,8 +21,9 @@ if __name__ == "__main__": parser.add_argument("--ch-jump", type=int) parser.add_argument("--ignore-flags", action="store_true") parser.add_argument("--autocorrelations-only", action="store_true", - help=("If this flag is set, all autocorrelations are imported " - "irrespective of whether they are flagged or not.")) + help=("If this flag is set, all autocorrelations are " + "imported irrespective of whether they are " + "flagged or not.")) parser.add_argument("ms") parser.add_argument( "polarization_mode", diff --git a/resolve/data/ms_import.py b/resolve/data/ms_import.py index aa7fe428..d089f299 100644 --- a/resolve/data/ms_import.py +++ b/resolve/data/ms_import.py @@ -275,19 +275,20 @@ def read_ms_i(name, data_column, freq, field, spectral_window, pol_indices, twgt = twgt[active_rows[start:stop]] twgt = twgt[:, active_channels] - tvis = t.getcol(data_column, startrow=start, nrow=stop - start)[ - ..., pol_indices - ] + # Vis + tvis = t.getcol(data_column, startrow=start, nrow=stop - start)[..., pol_indices] assert tvis.dtype == np.complex64 if not allrows: tvis = tvis[active_rows[start:stop]] tvis = tvis[:, active_channels] + # Flags tflags = _conditional_flags(t, start, stop, pol_indices, ignore_flags) if not allrows: tflags = tflags[active_rows[start:stop]] tflags = tflags[:, active_channels] + # Polarization summation assert twgt.ndim == tflags.ndim == 3 assert tflags.dtype == np.bool if not ignore_flags: @@ -316,15 +317,7 @@ def read_ms_i(name, data_column, freq, field, spectral_window, pol_indices, freq = freq[active_channels] my_asserteq(wgt.shape, vis.shape) - return ( - np.ascontiguousarray(uvw), - ant1, - ant2, - time, - np.ascontiguousarray(freq), - np.ascontiguousarray(vis), - np.ascontiguousarray(wgt), - ) + return {"uvw": uvw, "ant1": ant1, "ant2": ant2, "time": time, "freq": freq, "vis": vis, "wgt": wgt, "ptg": ptg} def ms_n_spectral_windows(ms): -- GitLab From 6008e8da4f9d0a2b300913ded29c86e5b12ff529 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 13 Aug 2021 11:32:47 +0200 Subject: [PATCH 16/89] How to deal with empty ms --- misc/ms2npz.py | 2 ++ resolve/data/ms_import.py | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/misc/ms2npz.py b/misc/ms2npz.py index 7c4fa716..0735f7bb 100755 --- a/misc/ms2npz.py +++ b/misc/ms2npz.py @@ -44,6 +44,8 @@ if __name__ == "__main__": args.ignore_flags ) for ifield, oo in enumerate(obs): + if oo is None: + continue if args.autocorrelations_only: oo = oo.restrict_to_autocorrelations() auto = "autocorrelationsonly" diff --git a/resolve/data/ms_import.py b/resolve/data/ms_import.py index d089f299..6ca2f20a 100644 --- a/resolve/data/ms_import.py +++ b/resolve/data/ms_import.py @@ -242,7 +242,7 @@ def read_ms_i(name, data_column, freq, field, spectral_window, pol_indices, start = stop nrealrows, nrealchan = np.sum(active_rows), np.sum(active_channels) if nrealrows == 0 or nrealchan == 0: - raise RuntimeError("Empty data set") + return None # Create output arrays if pol_summation: -- GitLab From 4f07aca134223e7afa621877251e8be57ad7197f Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 13 Aug 2021 11:33:05 +0200 Subject: [PATCH 17/89] Cosmetics --- resolve/data/ms_import.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/resolve/data/ms_import.py b/resolve/data/ms_import.py index 6ca2f20a..ec454729 100644 --- a/resolve/data/ms_import.py +++ b/resolve/data/ms_import.py @@ -210,9 +210,7 @@ def read_ms_i(name, data_column, freq, field, spectral_window, pol_indices, print("First pass:", f"{(start/nrow*100):.1f}%", end="\r") stop = min(nrow, start + step) tflags = _conditional_flags(t, start, stop, pol_indices, ignore_flags) - twgt = t.getcol(weightcol, startrow=start, nrow=stop - start)[ - ..., pol_indices - ] + twgt = t.getcol(weightcol, startrow=start, nrow=stop - start)[..., pol_indices] if channel_slice != slice(None): tchslcflags = np.ones_like(tflags) tchslcflags[:, channel_slice] = False -- GitLab From eb0d774055f9c3e7a5588b85b1cdf4ba5c6ff3d5 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 13 Aug 2021 11:33:28 +0200 Subject: [PATCH 18/89] Restructure and add pointing support --- misc/ms_summary.py | 5 ++ resolve/data/ms_import.py | 120 ++++++++++++++++++++++++-------------- 2 files changed, 81 insertions(+), 44 deletions(-) diff --git a/misc/ms_summary.py b/misc/ms_summary.py index 20f781c0..22550dd9 100755 --- a/misc/ms_summary.py +++ b/misc/ms_summary.py @@ -28,3 +28,8 @@ if __name__ == "__main__": print("NAME REFERNCE_DIR DELAY_DIR PHASE_DIR") for nn, rd, dd, pd in zip(name, refdir, deldir, phdir): print(nn, rd, dd, pd) + + with rve.ms_table(join(args.ms, "POINTING")) as t: + print(t.getcol("ANTENNA_ID")) + print(t.getcol("TIME")) + print(t.getcol("DIRECTION").shape) diff --git a/resolve/data/ms_import.py b/resolve/data/ms_import.py index ec454729..3ca3187d 100644 --- a/resolve/data/ms_import.py +++ b/resolve/data/ms_import.py @@ -84,8 +84,6 @@ def ms2observations(ms, data_column, with_calib_info, spectral_window, my_assert_isinstance(spectral_window, int) my_assert(spectral_window >= 0) my_assert(spectral_window < ms_n_spectral_windows(ms)) - with ms_table(join(ms, "SPECTRAL_WINDOW")) as t: - freq = t.getcol("CHAN_FREQ")[spectral_window] # Polarization with ms_table(join(ms, "POLARIZATION")) as t: @@ -110,11 +108,15 @@ def ms2observations(ms, data_column, with_calib_info, spectral_window, # Field with ms_table(join(ms, "FIELD")) as t: - equinox = t.coldesc("REFERENCE_DIR")["desc"]["keywords"]["MEASINFO"]["Ref"] - equinox = str(equinox)[1:] # FIXME Put proper support for equinox here - if equinox == "1950_VLA": - equinox = 1950 + # FIXME Eventually get rid of this and use auxiliary table + try: + equinox = t.coldesc("REFERENCE_DIR")["desc"]["keywords"]["MEASINFO"]["Ref"] + equinox = str(equinox)[1:] + if equinox == "1950_VLA": + equinox = 1950 + except KeyError: + equinox = 2000 dirs = [] for pc in t.getcol("REFERENCE_DIR"): my_asserteq(pc.shape, (1, 2)) @@ -125,37 +127,32 @@ def ms2observations(ms, data_column, with_calib_info, spectral_window, with ms_table(join(ms, "ANTENNA")) as t: keys = ["NAME", "STATION", "TYPE", "MOUNT", "POSITION", "OFFSET", "DISH_DIAMETER"] auxtables["ANTENNA"] = AuxiliaryTable({kk: t.getcol(kk) for kk in keys}) - with ms_table(join(ms, "FIELD")) as t: - keys = ["NAME", "CODE", "TIME", "NUM_POLY", "DELAY_DIR", "PHASE_DIR", "REFERENCE_DIR", - "SOURCE_ID"] - auxtables["FIELD"] = AuxiliaryTable({kk: t.getcol(kk) for kk in keys}) with ms_table(join(ms, "SPECTRAL_WINDOW")) as t: keys = ["NAME", "REF_FREQUENCY", "CHAN_FREQ", "CHAN_WIDTH", "MEAS_FREQ_REF", "EFFECTIVE_BW", "RESOLUTION", "TOTAL_BANDWIDTH", "NET_SIDEBAND", "IF_CONV_CHAIN", "FREQ_GROUP", "FREQ_GROUP_NAME"] - auxtables["SPECTRAL_WINDOW"] = AuxiliaryTable({kk: t.getcol(kk) for kk in keys}) + dct = {kk: t.getcol(kk, startrow=spectral_window, nrow=1) for kk in keys} + auxtables["SPECTRAL_WINDOW"] = AuxiliaryTable(dct) # FIXME Determine which observation is calibration observation - # FIXME Import name of source observations = [] for ifield, direction in enumerate(dirs): - uvw, ant1, ant2, time, freq_out, vis, wgt = read_ms_i( - ms, - data_column, - freq, - ifield, - spectral_window, - pol_ind, - pol_summation, - with_calib_info, - channel_slice, - ignore_flags, - ) - vis[wgt == 0] = 0.0 - vis = _ms2resolve_transpose(vis) - wgt = _ms2resolve_transpose(wgt) - antpos = AntennaPositions(uvw, ant1, ant2, time) - obs = Observation(antpos, vis, wgt, polobj, freq_out, direction, + with ms_table(join(ms, "FIELD")) as t: + keys = ["NAME", "CODE", "TIME", "NUM_POLY", "DELAY_DIR", "PHASE_DIR", "REFERENCE_DIR", + "SOURCE_ID"] + dct = {kk: t.getcol(kk, startrow=ifield, nrow=1) for kk in keys} + auxtables["FIELD"] = AuxiliaryTable(dct) + + mm = read_ms_i(ms, data_column, ifield, spectral_window, pol_ind, pol_summation, + with_calib_info, channel_slice, ignore_flags,) + if mm is None: + print(f"{ms}, field #{ifield} is empty or fully flagged") + observations.append(None) + continue + if mm["ptg"] is not None: + raise NotImplementedError + antpos = AntennaPositions(mm["uvw"], mm["ant1"], mm["ant2"], mm["time"]) + obs = Observation(antpos, mm["vis"], mm["wgt"], polobj, mm["freq"], direction, auxiliary_tables=auxtables) observations.append(obs) return observations @@ -178,12 +175,18 @@ def _determine_weighting(t): return fullwgt, weightcol -def read_ms_i(name, data_column, freq, field, spectral_window, pol_indices, - pol_summation, with_calib_info, channel_slice, ignore_flags): - freq = np.array(freq) - my_asserteq(freq.ndim, 1) +def read_ms_i(name, data_column, field, spectral_window, pol_indices, pol_summation, + with_calib_info, channel_slice, ignore_flags): + + # Freq + with ms_table(join(name, "SPECTRAL_WINDOW")) as t: + freq = t.getcol("CHAN_FREQ", startrow=spectral_window, nrow=1) + my_asserteq(freq.ndim, 2) + my_asserteq(freq.shape[0], 1) + freq = freq[0] my_assert(len(freq) > 0) nchan = len(freq) + assert pol_indices is None or isinstance(pol_indices, list) if pol_indices is None: pol_indices = slice(None) @@ -242,7 +245,11 @@ def read_ms_i(name, data_column, freq, field, spectral_window, pol_indices, if nrealrows == 0 or nrealchan == 0: return None - # Create output arrays + # Freq + freq = freq[active_channels] + + # Vis, wgt, (flags) + with ms_table(name) as t: if pol_summation: npol = 1 else: @@ -263,9 +270,8 @@ def read_ms_i(name, data_column, freq, field, spectral_window, pol_indices, if realstop > realstart: allrows = stop - start == realstop - realstart - twgt = t.getcol(weightcol, startrow=start, nrow=stop - start)[ - ..., pol_indices - ] + # Weights + twgt = t.getcol(weightcol, startrow=start, nrow=stop - start)[..., pol_indices] assert twgt.dtype == np.float32 if not fullwgt: twgt = np.repeat(twgt[:, None], nchan, axis=1) @@ -302,19 +308,45 @@ def read_ms_i(name, data_column, freq, field, spectral_window, pol_indices, wgt[realstart:realstop] = twgt start, realstart = stop, realstop - uvw = t.getcol("UVW")[active_rows] - if with_calib_info: + print("Selected:", 10 * " ") + print(f" shape: {vis.shape}") + print(f" flagged: {(1.0-np.sum(wgt!=0)/wgt.size)*100:.1f} %") + + # UVW + with ms_table(name) as t: + uvw = np.ascontiguousarray(t.getcol("UVW")[active_rows]) + + # Calibration info + if with_calib_info: + with ms_table(name) as t: ant1 = np.ascontiguousarray(t.getcol("ANTENNA1")[active_rows]) ant2 = np.ascontiguousarray(t.getcol("ANTENNA2")[active_rows]) time = np.ascontiguousarray(t.getcol("TIME")[active_rows]) + else: + ant1 = ant2 = time = None + + # Pointing + with ms_table(join(name, "POINTING")) as t: + if t.nrows() == 0: + ptg = None else: - ant1 = ant2 = time = None - print("Selected:", 10 * " ") - print(f" shape: {vis.shape}") - print(f" flagged: {(1.0-np.sum(wgt!=0)/wgt.size)*100:.1f} %") - freq = freq[active_channels] + ptg = np.empty((nrealrows, 1, 2), dtype=np.float64) + start, realstart = 0, 0 + while start < nrow: + print("Second pass:", f"{(start/nrow*100):.1f}%", end="\r") + stop = min(nrow, start + step) + realstop = realstart + np.sum(active_rows[start:stop]) + if realstop > realstart: + allrows = stop - start == realstop - realstart + tptg = t.getcol("DIRECTION", startrow=start, nrow=stop - start) + tptg = tptg[active_rows[start:stop]] + ptg[realstart:realstop] = tptg + start, realstart = stop, realstop my_asserteq(wgt.shape, vis.shape) + vis = np.ascontiguousarray(_ms2resolve_transpose(vis)) + wgt = np.ascontiguousarray(_ms2resolve_transpose(wgt)) + vis[wgt == 0] = 0.0 return {"uvw": uvw, "ant1": ant1, "ant2": ant2, "time": time, "freq": freq, "vis": vis, "wgt": wgt, "ptg": ptg} -- GitLab From 5819d35eb47206a27d80da6c37a8b909d6d5cef8 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 13 Aug 2021 11:57:20 +0200 Subject: [PATCH 19/89] NIFTy_7 -> NIFTy_8 --- misc/Dockerfile | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/misc/Dockerfile b/misc/Dockerfile index d7253334..82e729f4 100644 --- a/misc/Dockerfile +++ b/misc/Dockerfile @@ -2,8 +2,8 @@ FROM kernsuite/base:latest ENV DEBIAN_FRONTEND noninteractive RUN docker-apt-install dysco python3-casacore RUN apt-get update -qq && apt-get install -qq python3-pip python3-matplotlib git -RUN pip3 install scipy pybind11 git+https://gitlab.mpcdf.mpg.de/ift/nifty.git@NIFTy_7 -RUN pip3 install git+https://gitlab.mpcdf.mpg.de/mtr/ducc.git@ducc0 +RUN pip3 install scipy pybind11 git+https://gitlab.mpcdf.mpg.de/ift/nifty.git@NIFTy_8 +RUN pip3 install ducc0 COPY . /resolve RUN cd resolve && python3 setup.py install ENTRYPOINT ["python3", "resolve/misc/ms2npz.py"] -- GitLab From 5a541c08bfeeede8b6d8b75ab361277980b78bf1 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 26 Aug 2021 14:40:56 +0200 Subject: [PATCH 20/89] Data can be selected by frequency --- misc/ms2npz.py | 38 ++++++++++++++++++++++++++++++++++++-- misc/ms_summary.py | 14 +++++++++++--- resolve/constants.py | 27 +++++++++++++++++++++++++++ resolve/data/ms_import.py | 2 -- 4 files changed, 74 insertions(+), 7 deletions(-) diff --git a/misc/ms2npz.py b/misc/ms2npz.py index 0735f7bb..73505c28 100755 --- a/misc/ms2npz.py +++ b/misc/ms2npz.py @@ -4,9 +4,10 @@ # Author: Philipp Arras import argparse -from os.path import join, split, splitext from os import makedirs +from os.path import join, split, splitext +import numpy as np import resolve as rve if __name__ == "__main__": @@ -19,6 +20,8 @@ if __name__ == "__main__": parser.add_argument("--ch-begin", type=int) parser.add_argument("--ch-end", type=int) parser.add_argument("--ch-jump", type=int) + parser.add_argument("--freq-begin") + parser.add_argument("--freq-end") parser.add_argument("--ignore-flags", action="store_true") parser.add_argument("--autocorrelations-only", action="store_true", help=("If this flag is set, all autocorrelations are " @@ -30,6 +33,37 @@ if __name__ == "__main__": help="Can be 'stokesiavg', 'stokesi', 'all', or something like 'LL' or 'XY'.", ) args = parser.parse_args() + if args.freq_begin is not None and args.freq_end is not None: + assert args.ch_begin is None + assert args.ch_jump is None + assert args.ch_end is None + assert args.spectral_window == 0 + + f0 = rve.str2val(args.freq_begin) + f1 = rve.str2val(args.freq_end) + + # Determine spectral window + with rve.ms_table(join(args.ms, "SPECTRAL_WINDOW")) as t: + for ii in range(t.nrows()): + c = t.getcol("CHAN_FREQ", startrow=ii, nrow=1)[0] + fmin, fmax = min(c), max(c) + if fmin <= f0 <= fmax and fmin <= f1 <= fmax: + break + print(f"Load spectral window {ii}") + + # Determine channels + if np.all(np.diff(c) > 0): + begin = np.searchsorted(c, f0) + end = np.searchsorted(c, f1) + elif np.all(np.diff(c) < 0): + begin = c.size - np.searchsorted(c[::-1], f0) + end = c.size - np.searchsorted(c[::-1], f1) + else: + raise RuntimeError("Channels are not sorted") + channels = slice(begin, end) + else: + channels = slice(args.ch_begin, args.ch_end, args.ch_jump), + makedirs(args.output_folder, exist_ok=True) name = splitext(split(args.ms)[1])[0] nspec = rve.ms_n_spectral_windows(args.ms) @@ -40,7 +74,7 @@ if __name__ == "__main__": args.include_calibration_info, args.spectral_window, args.polarization_mode, - slice(args.ch_begin, args.ch_end, args.ch_jump), + channels, args.ignore_flags ) for ifield, oo in enumerate(obs): diff --git a/misc/ms_summary.py b/misc/ms_summary.py index 22550dd9..7a1734dd 100755 --- a/misc/ms_summary.py +++ b/misc/ms_summary.py @@ -15,9 +15,17 @@ if __name__ == "__main__": nspec = rve.ms_n_spectral_windows(args.ms) from casacore.tables import tablesummary tablesummary(args.ms) - print() - print(f"The data set has {nspec} spectral windows.") - print() + + with rve.ms_table(join(args.ms, "SPECTRAL_WINDOW")) as t: + for ii in range(nspec): + print(f"Spectral window #{ii}") + chans = t.getcol("CHAN_FREQ", startrow=ii, nrow=1) + + print(f"Shape: {chans.shape}") + print(f"f1-f0: {(chans[0][1]-chans[0][0])/1e6} MHz") + print("Frequencies (GHz)") + print(chans/1e9) + print() with rve.ms_table(join(args.ms, "FIELD")) as t: name = t.getcol("NAME") diff --git a/resolve/constants.py b/resolve/constants.py index ffd0999e..0a537533 100644 --- a/resolve/constants.py +++ b/resolve/constants.py @@ -38,3 +38,30 @@ def str2rad(s): if unit == kk: return float(s[:nn]) * c[kk] raise RuntimeError("Unit not understood") + + +def str2val(s): + """Convert string of number and unit to value. + + Support the following keys: p n mu m (nothing) k M G T + + Parameters + ---------- + s : str + TODO + + """ + c = { + "p": 1e-12, + "n": 1e-9, + "mu": 1e-6, + "m": 1e-3, + "k": 1e3, + "M": 1e6, + "G": 1e9, + "T": 1e12 + } + keys = set(c.keys()) + if s[-1] in keys: + return float(s[:-1]) * c[s[-1]] + return float(s) diff --git a/resolve/data/ms_import.py b/resolve/data/ms_import.py index 3ca3187d..1dd4310b 100644 --- a/resolve/data/ms_import.py +++ b/resolve/data/ms_import.py @@ -68,8 +68,6 @@ def ms2observations(ms, data_column, with_calib_info, spectral_window, raise RuntimeError if ms == ".": ms = os.getcwd() - if splitext(ms)[1].lower() != ".ms": - raise RuntimeError if isinstance(polarizations, str): polarizations = [polarizations] if ( -- GitLab From 57dae67b922a9c52968ea6157789db8c4bf30656 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 26 Aug 2021 18:00:40 +0200 Subject: [PATCH 21/89] FIXME --- resolve/data/ms_import.py | 1 + 1 file changed, 1 insertion(+) diff --git a/resolve/data/ms_import.py b/resolve/data/ms_import.py index 1dd4310b..073f76a2 100644 --- a/resolve/data/ms_import.py +++ b/resolve/data/ms_import.py @@ -139,6 +139,7 @@ def ms2observations(ms, data_column, with_calib_info, spectral_window, keys = ["NAME", "CODE", "TIME", "NUM_POLY", "DELAY_DIR", "PHASE_DIR", "REFERENCE_DIR", "SOURCE_ID"] dct = {kk: t.getcol(kk, startrow=ifield, nrow=1) for kk in keys} + # FIXME Somehow not the correct field is loaded here auxtables["FIELD"] = AuxiliaryTable(dct) mm = read_ms_i(ms, data_column, ifield, spectral_window, pol_ind, pol_summation, -- GitLab From 6ad51f87233d2f8f6b6b5041849cf1c3150f5a1b Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 27 Aug 2021 11:32:50 +0200 Subject: [PATCH 22/89] Correctly read in polarization index --- resolve/data/ms_import.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/resolve/data/ms_import.py b/resolve/data/ms_import.py index 073f76a2..67e502cb 100644 --- a/resolve/data/ms_import.py +++ b/resolve/data/ms_import.py @@ -83,12 +83,15 @@ def ms2observations(ms, data_column, with_calib_info, spectral_window, my_assert(spectral_window >= 0) my_assert(spectral_window < ms_n_spectral_windows(ms)) + with ms_table(join(ms, "DATA_DESCRIPTION")) as t: + polid = t.getcol("POLARIZATION_ID")[spectral_window] + # Polarization with ms_table(join(ms, "POLARIZATION")) as t: - pol = t.getcol("CORR_TYPE") + pol = t.getcol("CORR_TYPE", startrow=polid, nrow=1) my_asserteq(pol.ndim, 2) my_asserteq(pol.shape[0], 1) - pol = list(pol[0]) # Not clear what the first dimension is used for + pol = list(pol[0]) polobj = Polarization(pol) if polarizations[0] == "stokesi": polarizations = ["LL", "RR"] if polobj.circular() else ["XX", "YY"] -- GitLab From 443373372c770bd4cfc628a68785d1e5a7bd1cd6 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 27 Aug 2021 11:41:44 +0200 Subject: [PATCH 23/89] Fixup --- resolve/plotter.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/resolve/plotter.py b/resolve/plotter.py index 598a7081..9d4de344 100644 --- a/resolve/plotter.py +++ b/resolve/plotter.py @@ -203,7 +203,8 @@ def _plot_nifty(state, op, kwargs, fname): sc.add(op.force(ss)) p = ift.Plot() p.add(sc.mean, **kwargs) - p.add(sc.var.sqrt() / sc.mean) + if len(state) > 1: + p.add(sc.var.sqrt() / sc.mean) p.output(nx=2, ny=1, xsize=2 * UNIT, ysize=UNIT, name=fname) else: pos = state if isinstance(state, ift.MultiField) else state.mean -- GitLab From f7c0476171116d2c0b05bb1d141536741bbd885d Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Mon, 30 Aug 2021 17:12:25 +0200 Subject: [PATCH 24/89] Cleanup ms import and use auxiliary tables more --- misc/ms2npz.py | 2 +- resolve/data/auxiliary_table.py | 12 ++- resolve/data/ms_import.py | 161 +++++++++++++++++--------------- resolve/data/observation.py | 39 ++++---- resolve/util.py | 8 +- test/test_general.py | 9 +- 6 files changed, 127 insertions(+), 104 deletions(-) diff --git a/misc/ms2npz.py b/misc/ms2npz.py index 73505c28..7d8f427c 100755 --- a/misc/ms2npz.py +++ b/misc/ms2npz.py @@ -62,7 +62,7 @@ if __name__ == "__main__": raise RuntimeError("Channels are not sorted") channels = slice(begin, end) else: - channels = slice(args.ch_begin, args.ch_end, args.ch_jump), + channels = slice(args.ch_begin, args.ch_end, args.ch_jump) makedirs(args.output_folder, exist_ok=True) name = splitext(split(args.ms)[1])[0] diff --git a/resolve/data/auxiliary_table.py b/resolve/data/auxiliary_table.py index 58365b8c..a2afcb71 100644 --- a/resolve/data/auxiliary_table.py +++ b/resolve/data/auxiliary_table.py @@ -13,7 +13,8 @@ class AuxiliaryTable: nrow = None for kk, lst in inp_dict.items(): my_assert_isinstance(kk, str) - my_assert_isinstance(lst, (list, np.ndarray)) + if not isinstance(lst, (list, np.ndarray)): + raise RuntimeError(f"{kk} neither list nor np.ndarray") if nrow is None: nrow = len(lst) my_asserteq(nrow, len(lst)) @@ -35,6 +36,15 @@ class AuxiliaryTable: s.append(f" {kk:<20} {len(lst):>10} {str(type(lst[0])):>15}") return "\n".join(s) + def __getitem__(self, key): + return self._dct[key] + + def __contains__(self, key): + return key in self._dct + + def row(self, i): + return AuxiliaryTable({kk: vv[i:i+1] for kk, vv in self._dct.items()}) + def to_list(self): return [list(self._dct.keys())] + [ss for ss in self._dct.values()] diff --git a/resolve/data/ms_import.py b/resolve/data/ms_import.py index 67e502cb..e3214afd 100644 --- a/resolve/data/ms_import.py +++ b/resolve/data/ms_import.py @@ -19,9 +19,16 @@ def ms_table(path): return table(path, readonly=True, ack=False) +def _pol_id(ms_path, spectral_window): + """Return id for indexing polarization table for a given spectral window.""" + with ms_table(join(ms_path, "DATA_DESCRIPTION")) as t: + polid = t.getcol("POLARIZATION_ID")[spectral_window] + return polid + + + def ms2observations(ms, data_column, with_calib_info, spectral_window, - polarizations="all", channel_slice=slice(None), - ignore_flags=False): + polarizations="all", channels=slice(None), ignore_flags=False): """Read and convert a given measurement set into an array of :class:`Observation` If WEIGHT_SPECTRUM is available this column is used for weighting. @@ -38,15 +45,15 @@ def ms2observations(ms, data_column, with_calib_info, spectral_window, Reads in all information necessary for calibration, if True. If only imaging shall be performed, this can be set to False. spectral_window : int - Index of spectral window which shall be imported. + Index of spectral window that shall be imported. polarizations "all": All polarizations are imported. "stokesi": Only LL/RR or XX/YY polarizations are imported. "stokesiavg": Only LL/RR or XX/YY polarizations are imported and averaged on the fly. List of strings: Strings can be "XX", "YY", "XY", "YX", "LL", "LR", "RL", "RR". The respective polarizations are loaded. - channel_slice : slice - Slice of selected channels. Default: select all channels - FIXME Select channels by indices + channels : slice or list + Select channels. Can be either a slice object or an index list. Default: + select all channels ignore_flags : bool If True, the whole measurement set is imported irrespective of the flags. Default is false. @@ -58,95 +65,63 @@ def ms2observations(ms, data_column, with_calib_info, spectral_window, Note ---- - We cannot import multiple spectral windows into one Observation instance + Multiple spectral windows are not imported into one Observation instance because it is not guaranteed by the measurement set data structure that all baselines are present in all spectral windows. """ + # Input checks if ms[-1] == "/": ms = ms[:-1] if not isdir(ms): raise RuntimeError if ms == ".": ms = os.getcwd() - if isinstance(polarizations, str): - polarizations = [polarizations] - if ( - "stokesiavg" in polarizations - or "stokesi" in polarizations - or "all" in polarizations - ) and len(polarizations) > 1: - raise ValueError - - # Spectral windows - my_assert_isinstance(channel_slice, slice) + if isinstance(polarizations, list): + for ll in polarizations: + my_assert(ll in ["XX", "YY", "XY", "YX", "LL", "LR", "RL", "RR"]) + my_asserteq(len(set(polarizations)), len(polarizations)) + else: + my_assert(polarizations in ["stokesi", "stokesiavg", "all"]) + my_assert_isinstance(channels, (slice, list)) my_assert_isinstance(spectral_window, int) my_assert(spectral_window >= 0) my_assert(spectral_window < ms_n_spectral_windows(ms)) - - with ms_table(join(ms, "DATA_DESCRIPTION")) as t: - polid = t.getcol("POLARIZATION_ID")[spectral_window] + # /Input checks # Polarization with ms_table(join(ms, "POLARIZATION")) as t: - pol = t.getcol("CORR_TYPE", startrow=polid, nrow=1) - my_asserteq(pol.ndim, 2) - my_asserteq(pol.shape[0], 1) - pol = list(pol[0]) + pol = t.getcol("CORR_TYPE", startrow=_pol_id(ms, spectral_window), nrow=1)[0] polobj = Polarization(pol) - if polarizations[0] == "stokesi": - polarizations = ["LL", "RR"] if polobj.circular() else ["XX", "YY"] - if polarizations[0] == "stokesiavg": + if polarizations == "stokesiavg": pol_ind = polobj.stokes_i_indices() polobj = Polarization.trivial() pol_summation = True - elif polarizations[0] == "all": + elif polarizations == "all": pol_ind = None pol_summation = False else: + if polarizations == "stokesi": + polarizations = ["LL", "RR"] if polobj.circular() else ["XX", "YY"] polobj = polobj.restrict_by_name(polarizations) pol_ind = [polobj.to_str_list().index(ii) for ii in polarizations] pol_summation = False - # Field - with ms_table(join(ms, "FIELD")) as t: - # FIXME Put proper support for equinox here - # FIXME Eventually get rid of this and use auxiliary table - try: - equinox = t.coldesc("REFERENCE_DIR")["desc"]["keywords"]["MEASINFO"]["Ref"] - equinox = str(equinox)[1:] - if equinox == "1950_VLA": - equinox = 1950 - except KeyError: - equinox = 2000 - dirs = [] - for pc in t.getcol("REFERENCE_DIR"): - my_asserteq(pc.shape, (1, 2)) - dirs.append(Direction(pc[0], equinox)) - dirs = tuple(dirs) - - auxtables = {} - with ms_table(join(ms, "ANTENNA")) as t: - keys = ["NAME", "STATION", "TYPE", "MOUNT", "POSITION", "OFFSET", "DISH_DIAMETER"] - auxtables["ANTENNA"] = AuxiliaryTable({kk: t.getcol(kk) for kk in keys}) - with ms_table(join(ms, "SPECTRAL_WINDOW")) as t: - keys = ["NAME", "REF_FREQUENCY", "CHAN_FREQ", "CHAN_WIDTH", "MEAS_FREQ_REF", "EFFECTIVE_BW", - "RESOLUTION", "TOTAL_BANDWIDTH", "NET_SIDEBAND", "IF_CONV_CHAIN", "FREQ_GROUP", - "FREQ_GROUP_NAME"] - dct = {kk: t.getcol(kk, startrow=spectral_window, nrow=1) for kk in keys} - auxtables["SPECTRAL_WINDOW"] = AuxiliaryTable(dct) - - # FIXME Determine which observation is calibration observation observations = [] - for ifield, direction in enumerate(dirs): - with ms_table(join(ms, "FIELD")) as t: - keys = ["NAME", "CODE", "TIME", "NUM_POLY", "DELAY_DIR", "PHASE_DIR", "REFERENCE_DIR", - "SOURCE_ID"] - dct = {kk: t.getcol(kk, startrow=ifield, nrow=1) for kk in keys} - # FIXME Somehow not the correct field is loaded here - auxtables["FIELD"] = AuxiliaryTable(dct) + for ifield in range(ms_n_fields(ms)): + auxtables = {} + auxtables["ANTENNA"] = _import_aux_table(ms, "ANTENNA") + auxtables["STATE"] = _import_aux_table(ms, "STATE") + auxtables["SPECTRAL_WINDOW"] = _import_aux_table(ms, "SPECTRAL_WINDOW", row=spectral_window, skip=["ASSOC_NATURE"]) + sf = _source_and_field_table(ms, spectral_window, ifield) + if sf is None: + print(f"Field {ifield} cannot be found in SOURCE table") + observations.append(None) + continue + auxtables = {**auxtables, **sf} + print(f"Work on Field {ifield}: {auxtables['SOURCE']['NAME'][0]}") mm = read_ms_i(ms, data_column, ifield, spectral_window, pol_ind, pol_summation, - with_calib_info, channel_slice, ignore_flags,) + with_calib_info, channels, ignore_flags) if mm is None: print(f"{ms}, field #{ifield} is empty or fully flagged") observations.append(None) @@ -154,7 +129,7 @@ def ms2observations(ms, data_column, with_calib_info, spectral_window, if mm["ptg"] is not None: raise NotImplementedError antpos = AntennaPositions(mm["uvw"], mm["ant1"], mm["ant2"], mm["time"]) - obs = Observation(antpos, mm["vis"], mm["wgt"], polobj, mm["freq"], direction, + obs = Observation(antpos, mm["vis"], mm["wgt"], polobj, mm["freq"], auxiliary_tables=auxtables) observations.append(obs) return observations @@ -178,14 +153,12 @@ def _determine_weighting(t): def read_ms_i(name, data_column, field, spectral_window, pol_indices, pol_summation, - with_calib_info, channel_slice, ignore_flags): + with_calib_info, channels, ignore_flags): # Freq with ms_table(join(name, "SPECTRAL_WINDOW")) as t: - freq = t.getcol("CHAN_FREQ", startrow=spectral_window, nrow=1) - my_asserteq(freq.ndim, 2) - my_asserteq(freq.shape[0], 1) - freq = freq[0] + freq = t.getcol("CHAN_FREQ", startrow=spectral_window, nrow=1)[0] + my_asserteq(freq.ndim, 1) my_assert(len(freq) > 0) nchan = len(freq) @@ -216,10 +189,11 @@ def read_ms_i(name, data_column, field, spectral_window, pol_indices, pol_summat stop = min(nrow, start + step) tflags = _conditional_flags(t, start, stop, pol_indices, ignore_flags) twgt = t.getcol(weightcol, startrow=start, nrow=stop - start)[..., pol_indices] - if channel_slice != slice(None): - tchslcflags = np.ones_like(tflags) - tchslcflags[:, channel_slice] = False - tflags = np.logical_or(tflags, tchslcflags) + + tchslcflags = np.ones_like(tflags) + tchslcflags[:, channels] = False + tflags = np.logical_or(tflags, tchslcflags) + if not fullwgt: twgt = np.repeat(twgt[:, None], nchan, axis=1) my_asserteq(twgt.ndim, 3) @@ -358,8 +332,43 @@ def ms_n_spectral_windows(ms): return n_spectral_windows +def ms_n_fields(ms): + with ms_table(join(ms, "FIELD")) as t: + n = t.nrows() + return n + + def _conditional_flags(table, start, stop, pol_indices, ignore): tflags = table.getcol("FLAG", startrow=start, nrow=stop - start)[..., pol_indices] if ignore: tflags = np.zeros_like(tflags) return tflags + + +def _import_aux_table(ms, table_name, row=None, skip=[]): + with ms_table(join(ms, table_name)) as t: + keys = filter(lambda s: s not in skip, t.colnames()) + if row is None: + dct = {kk: t.getcol(kk) for kk in keys} + else: + dct = {kk: t.getcol(kk, startrow=row, nrow=1) for kk in keys} + aux_table = AuxiliaryTable(dct) + return aux_table + + +def _source_and_field_table(ms, spectral_window, ifield): + source_table = _import_aux_table(ms, "SOURCE", skip=["POSITION", "TRANSITION", "REST_FREQUENCY", + "SYSVEL", "SOURCE_MODEL", "PULSAR_ID"]) + field_table = _import_aux_table(ms, "FIELD", row=ifield) + source_id = field_table["SOURCE_ID"][0] + ind = np.where(np.logical_and(source_table["SOURCE_ID"] == source_id, + np.logical_or(source_table["SPECTRAL_WINDOW_ID"] == spectral_window, + source_table["SPECTRAL_WINDOW_ID"] == -1)))[0] + if len(ind) == 0: + return None + elif len(ind) == 1: + ind = ind[0] + else: + raise RuntimeError + + return {"SOURCE": source_table.row(ind), "FIELD": field_table} diff --git a/resolve/data/observation.py b/resolve/data/observation.py index 8e035487..73407d08 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -267,8 +267,6 @@ class Observation(BaseObservation): Polarization information of the data set. freq : numpy.ndarray Contains the measured frequencies. Shape (n_channels) - direction : Direction - Direction information of the data set. auxiliary_tables : dict Dictionary of auxiliary tables. Default: None. @@ -277,9 +275,8 @@ class Observation(BaseObservation): vis and weight must have the same shape. """ - def __init__(self, antenna_positions, vis, weight, polarization, freq, direction, auxiliary_tables=None): + def __init__(self, antenna_positions, vis, weight, polarization, freq, auxiliary_tables=None): nrows = len(antenna_positions) - my_assert_isinstance(direction, Direction) my_assert_isinstance(polarization, Polarization) my_assert_isinstance(antenna_positions, AntennaPositions) my_asserteq(weight.shape, vis.shape) @@ -305,7 +302,6 @@ class Observation(BaseObservation): self._weight = weight self._polarization = polarization self._freq = freq - self._direction = direction if auxiliary_tables is not None: my_assert_isinstance(auxiliary_tables, dict) @@ -314,8 +310,8 @@ class Observation(BaseObservation): my_assert_isinstance(kk, str) self._auxiliary_tables = auxiliary_tables - self._eq_attributes = ("_direction", "_polarization", "_freq", - "_antpos", "_vis", "_weight", "_auxiliary_tables") + self._eq_attributes = ("_polarization", "_freq", "_antpos", "_vis", "_weight", + "_auxiliary_tables") @onlymaster def save(self, file_name, compress): @@ -324,7 +320,6 @@ class Observation(BaseObservation): weight=self._weight, freq=self._freq, polarization=self._polarization.to_list(), - direction=self._direction.to_list(), ) for ii, vv in enumerate(self._antpos.to_list()): if vv is None: @@ -363,7 +358,6 @@ class Observation(BaseObservation): # /Load auxtables pol = Polarization.from_list(dct["polarization"]) - direction = Direction.from_list(dct["direction"]) slc = slice(None) if lo_hi_index is None else slice(*lo_hi_index) # FIXME Put barrier here that makes sure that only one full Observation is loaded at a time return Observation( @@ -372,7 +366,6 @@ class Observation(BaseObservation): dct["weight"][..., slc], pol, dct["freq"][slc], - direction, auxtables ) @@ -382,7 +375,7 @@ class Observation(BaseObservation): vis = self._vis.copy() vis[self.flags.val] = np.nan return Observation(self._antpos, vis, self._weight, self._polarization, self._freq, - self._direction) + self._auxiliary_tables) @staticmethod def load_mf(file_name, n_imaging_bands, comm=None): @@ -413,13 +406,14 @@ class Observation(BaseObservation): return obs_list, nu0 def __getitem__(self, slc): + # FIXME Do I need to change something in self._auxiliary_tables? return Observation( self._antpos[slc], self._vis[:, slc], self._weight[:, slc], self._polarization, self._freq, - self._direction, + self._auxiliary_tables, ) def get_freqs(self, frequency_list): @@ -435,16 +429,18 @@ class Observation(BaseObservation): return self.get_freqs_by_slice(mask) def get_freqs_by_slice(self, slc): + # FIXME Do I need to change something in self._auxiliary_tables? return Observation( self._antpos, self._vis[..., slc], self._weight[..., slc], self._polarization, self._freq[slc], - self._direction, + self._auxiliary_tables, ) def average_stokesi(self): + # FIXME Do I need to change something in self._auxiliary_tables? my_asserteq(self._vis.shape[0], 2) my_asserteq(self._polarization.restrict_to_stokes_i(), self._polarization) vis = np.sum(self._weight * self._vis, axis=0)[None] @@ -453,22 +449,24 @@ class Observation(BaseObservation): vis /= wgt + np.ones_like(wgt) * invmask vis[invmask] = 0.0 return Observation( - self._antpos, vis, wgt, Polarization.trivial(), self._freq, self._direction + self._antpos, vis, wgt, Polarization.trivial(), self._freq, self._auxiliary_tables ) def restrict_to_stokesi(self): + # FIXME Do I need to change something in self._auxiliary_tables? my_asserteq(self._vis.shape[0], 4) ind = self._polarization.stokes_i_indices() vis = self._vis[ind] wgt = self._weight[ind] pol = self._polarization.restrict_to_stokes_i() - return Observation(self._antpos, vis, wgt, pol, self._freq, self._direction) + return Observation(self._antpos, vis, wgt, pol, self._freq, self._auxiliary_tables) def restrict_to_autocorrelations(self): slc = self._antpos.ant1 == self._antpos.ant2 return self[slc] def move_time(self, t0): + # FIXME Do I need to change something in self._auxiliary_tables? antpos = self._antpos.move_time(t0) return Observation( antpos, @@ -476,7 +474,7 @@ class Observation(BaseObservation): self._weight, self._polarization, self._freq, - self._direction, + self._auxiliary_tables ) @property @@ -487,6 +485,15 @@ class Observation(BaseObservation): def antenna_positions(self): return self._antpos + @property + def direction(self): + if self._auxiliary_tables is not None and "FIELD" in self._auxiliary_tables: + equinox = 2000 # FIXME Figure out how to extract this from a measurement set + refdir = self._auxiliary_tables["FIELD"]["REFERENCE_DIR"][0] + my_asserteq(refdir.shape[0] == 0) + return Direction(refdir[0], equinox) + return None + def auxiliary_table(self, name): return self._auxiliary_tables[name] diff --git a/resolve/util.py b/resolve/util.py index b8eda2b3..8093ef2e 100644 --- a/resolve/util.py +++ b/resolve/util.py @@ -60,10 +60,10 @@ def _deep_equal(a, b): def _compare_simple_or_array(a, b): equal = a == b - try: - return bool(equal) - except ValueError: - return equal.all() + if isinstance(equal, np.ndarray): + return np.all(equal) + assert isinstance(equal, bool) + return equal def _types_equal(a, b): diff --git a/test/test_general.py b/test/test_general.py index e5e2f997..13845c50 100644 --- a/test/test_general.py +++ b/test/test_general.py @@ -42,10 +42,7 @@ freqdomain = rve.IRGSpace(np.linspace(1000, 1050, num=10)) FACETS = [(1, 1), (2, 2), (2, 1), (1, 4)] -@pmp( - "ms", - ("CYG-ALL-2052-2MHZ", "CYG-D-6680-64CH-10S", "AM754_A030124_flagged", "1052735056"), -) +@pmp("ms", ("CYG-ALL-2052-2MHZ", "CYG-D-6680-64CH-10S", "AM754_A030124_flagged")) @pmp("with_calib_info", (False, True)) @pmp("compress", (False, True)) def test_save_and_load_observation(ms, with_calib_info, compress): @@ -58,11 +55,11 @@ def test_save_and_load_observation(ms, with_calib_info, compress): assert ob == ob1 -@pmp("slc", [slice(3), slice(14, 15), slice(None, None, None), slice(None, None, 5)]) +@pmp("slc", [slice(3), slice(14, 15), [14, 15], slice(None, None, None), slice(None, None, 5)]) def test_sliced_readin(slc): ms = f"{direc}CYG-D-6680-64CH-10S.ms" obs0 = rve.ms2observations(ms, "DATA", False, 0)[0] - obs = rve.ms2observations(ms, "DATA", False, 0, channel_slice=slc)[0] + obs = rve.ms2observations(ms, "DATA", False, 0, channels=slc)[0] ch0 = obs0.weight.val[..., slc] ch = obs.weight.val assert ch0.ndim == 3 -- GitLab From 443fcff27ef6a2174b8a190d602f8e54995d0d60 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Tue, 31 Aug 2021 13:01:35 +0200 Subject: [PATCH 25/89] ms_import.py: Refactoring --- resolve/data/ms_import.py | 174 ++++++++++++++++++++++---------------- 1 file changed, 99 insertions(+), 75 deletions(-) diff --git a/resolve/data/ms_import.py b/resolve/data/ms_import.py index e3214afd..c764a9e7 100644 --- a/resolve/data/ms_import.py +++ b/resolve/data/ms_import.py @@ -2,14 +2,13 @@ # Copyright(C) 2019-2021 Max-Planck-Society import os -from os.path import isdir, join, splitext +from os.path import isdir, join import numpy as np from ..util import my_assert, my_assert_isinstance, my_asserteq from .antenna_positions import AntennaPositions from .auxiliary_table import AuxiliaryTable -from .direction import Direction from .observation import Observation from .polarization import Polarization @@ -26,7 +25,6 @@ def _pol_id(ms_path, spectral_window): return polid - def ms2observations(ms, data_column, with_calib_info, spectral_window, polarizations="all", channels=slice(None), ignore_flags=False): """Read and convert a given measurement set into an array of :class:`Observation` @@ -140,91 +138,44 @@ def _ms2resolve_transpose(arr): return np.ascontiguousarray(np.transpose(arr, (2, 0, 1))) -def _determine_weighting(t): - fullwgt = False - weightcol = "WEIGHT" - try: - t.getcol("WEIGHT_SPECTRUM", startrow=0, nrow=1) - weightcol = "WEIGHT_SPECTRUM" - fullwgt = True - except RuntimeError: - pass +def _determine_weighting(ms): + with ms_table(ms) as t: + if "WEIGHT_SPECTRUM" in t.colnames(): + weightcol = "WEIGHT_SPECTRUM" + fullwgt = True + else: + weightcol = "WEIGHT" + fullwgt = False return fullwgt, weightcol def read_ms_i(name, data_column, field, spectral_window, pol_indices, pol_summation, with_calib_info, channels, ignore_flags): - - # Freq - with ms_table(join(name, "SPECTRAL_WINDOW")) as t: - freq = t.getcol("CHAN_FREQ", startrow=spectral_window, nrow=1)[0] - my_asserteq(freq.ndim, 1) - my_assert(len(freq) > 0) - nchan = len(freq) - assert pol_indices is None or isinstance(pol_indices, list) if pol_indices is None: pol_indices = slice(None) if pol_summation: my_asserteq(len(pol_indices), 2) + # Check if data column is available and get shape with ms_table(name) as t: - # FIXME Get rid of fullwgt - fullwgt, weightcol = _determine_weighting(t) - nrow = t.nrows() - nmspol = t.getcol("FLAG", startrow=0, nrow=1).shape[2] - print("Measurement set visibilities:") - print(f" shape: ({nrow}, {nchan}, {nmspol})") - active_rows = np.ones(nrow, dtype=np.bool) - active_channels = np.zeros(nchan, dtype=np.bool) - step = max(1, nrow // 100) # how many rows to read in every step - - # Check if data column is available - t.getcol(data_column, startrow=0, nrow=10) - - # Determine which subset of rows/channels we need to input - start = 0 - while start < nrow: - print("First pass:", f"{(start/nrow*100):.1f}%", end="\r") - stop = min(nrow, start + step) - tflags = _conditional_flags(t, start, stop, pol_indices, ignore_flags) - twgt = t.getcol(weightcol, startrow=start, nrow=stop - start)[..., pol_indices] - - tchslcflags = np.ones_like(tflags) - tchslcflags[:, channels] = False - tflags = np.logical_or(tflags, tchslcflags) - - if not fullwgt: - twgt = np.repeat(twgt[:, None], nchan, axis=1) - my_asserteq(twgt.ndim, 3) - npol = tflags.shape[2] - if pol_summation: - tflags = np.any(tflags.astype(np.bool), axis=2)[..., None] - twgt = np.sum(twgt, axis=2)[..., None] - tflags[twgt == 0] = True - - # Select field and spectral window - tfieldid = t.getcol("FIELD_ID", startrow=start, nrow=stop - start) - tflags[tfieldid != field] = True - tspw = t.getcol("DATA_DESC_ID", startrow=start, nrow=stop - start) - tflags[tspw != spectral_window] = True - - # Inactive if all polarizations are flagged - assert tflags.ndim == 3 - tflags = np.all(tflags, axis=2) - active_rows[start:stop] = np.invert(np.all(tflags, axis=1)) - active_channels = np.logical_or( - active_channels, np.invert(np.all(tflags, axis=0)) - ) - start = stop - nrealrows, nrealchan = np.sum(active_rows), np.sum(active_channels) - if nrealrows == 0 or nrealchan == 0: - return None + nmspol = t.getcol(data_column, startrow=0, nrow=3).shape[2] + nrow = t.nrow() + print("Measurement set visibilities:") + print(f" shape: ({nrow}, {_ms_nchannels(name, spectral_window)}, {nmspol})") + + active_rows, active_channels = _first_pass(name, field, spectral_window, channels, pol_indices, + pol_summation, ignore_flags) + nrealrows, nrealchan = np.sum(active_rows), np.sum(active_channels) + if nrealrows == 0 or nrealchan == 0: + return None # Freq - freq = freq[active_channels] + freq = _ms_channels(name, spectral_window)[active_channels] # Vis, wgt, (flags) + fullwgt, weightcol = _determine_weighting(name) + nchan = _ms_nchannels(name, spectral_window) with ms_table(name) as t: if pol_summation: npol = 1 @@ -241,7 +192,7 @@ def read_ms_i(name, data_column, field, spectral_window, pol_indices, pol_summat start, realstart = 0, 0 while start < nrow: print("Second pass:", f"{(start/nrow*100):.1f}%", end="\r") - stop = min(nrow, start + step) + stop = _ms_stop(start, nrow) realstop = realstart + np.sum(active_rows[start:stop]) if realstop > realstart: allrows = stop - start == realstop - realstart @@ -310,7 +261,7 @@ def read_ms_i(name, data_column, field, spectral_window, pol_indices, pol_summat start, realstart = 0, 0 while start < nrow: print("Second pass:", f"{(start/nrow*100):.1f}%", end="\r") - stop = min(nrow, start + step) + stop = _ms_stop(start, nrow) realstop = realstart + np.sum(active_rows[start:stop]) if realstop > realstart: allrows = stop - start == realstop - realstart @@ -323,7 +274,80 @@ def read_ms_i(name, data_column, field, spectral_window, pol_indices, pol_summat vis = np.ascontiguousarray(_ms2resolve_transpose(vis)) wgt = np.ascontiguousarray(_ms2resolve_transpose(wgt)) vis[wgt == 0] = 0.0 - return {"uvw": uvw, "ant1": ant1, "ant2": ant2, "time": time, "freq": freq, "vis": vis, "wgt": wgt, "ptg": ptg} + return {"uvw": uvw, "ant1": ant1, "ant2": ant2, "time": time, "freq": freq, "vis": vis, + "wgt": wgt, "ptg": ptg} + + +def _first_pass(ms, field, spectral_window, channels, pol_indices, pol_summation, ignore_flags): + """Go through measurement set and determine which rows and which channels are active for a given + field and a given spectral window. + """ + fullwgt, weightcol = _determine_weighting(ms) + nchan = _ms_nchannels(ms, spectral_window) + with ms_table(ms) as t: + nrow = t.nrows() + active_rows = np.ones(nrow, dtype=np.bool) + active_channels = np.zeros(nchan, dtype=np.bool) + + # Determine which subset of rows/channels we need to input + start = 0 + while start < nrow: + print("First pass:", f"{(start/nrow*100):.1f}%", end="\r") + stop = _ms_stop(start, nrow) + tflags = _conditional_flags(t, start, stop, pol_indices, ignore_flags) + twgt = t.getcol(weightcol, startrow=start, nrow=stop - start)[..., pol_indices] + + tchslcflags = np.ones_like(tflags) + tchslcflags[:, channels] = False + tflags = np.logical_or(tflags, tchslcflags) + + if not fullwgt: + twgt = np.repeat(twgt[:, None], nchan, axis=1) + my_asserteq(twgt.ndim, 3) + if pol_summation: + tflags = np.any(tflags.astype(np.bool), axis=2)[..., None] + twgt = np.sum(twgt, axis=2)[..., None] + tflags[twgt == 0] = True + + # Select field and spectral window + tfieldid = t.getcol("FIELD_ID", startrow=start, nrow=stop - start) + tflags[tfieldid != field] = True + tspw = t.getcol("DATA_DESC_ID", startrow=start, nrow=stop - start) + tflags[tspw != spectral_window] = True + + # Inactive if all polarizations are flagged + assert tflags.ndim == 3 + tflags = np.all(tflags, axis=2) + active_rows[start:stop] = np.invert(np.all(tflags, axis=1)) + active_channels = np.logical_or( + active_channels, np.invert(np.all(tflags, axis=0)) + ) + start = stop + return active_rows, active_channels + + +def _ms_stop(start, nrow): + """Compute sensible step size for going through the rows of a measurement set and return stop + index. + """ + step = max(1, nrow // 100) + return min(nrow, start + step) + + +def _ms_nchannels(ms, spectral_window): + """Return number of channels in a given spectral window of a measurement set. + """ + return len(_ms_channels(ms, spectral_window)) + + +def _ms_channels(ms, spectral_window): + """Return frequencies of channels in a given spectral window of a measurement set. + """ + with ms_table(join(ms, "SPECTRAL_WINDOW")) as t: + freq = t.getcol("CHAN_FREQ", startrow=spectral_window, nrow=1)[0] + my_asserteq(freq.ndim, 1) + my_assert(len(freq) > 0) + return freq def ms_n_spectral_windows(ms): -- GitLab From c0ed2537c61357c85978a160d4f891a050af812c Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Tue, 31 Aug 2021 13:04:29 +0200 Subject: [PATCH 26/89] fixup! ms_import.py: Refactoring --- resolve/data/ms_import.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/resolve/data/ms_import.py b/resolve/data/ms_import.py index c764a9e7..0678a31c 100644 --- a/resolve/data/ms_import.py +++ b/resolve/data/ms_import.py @@ -160,7 +160,7 @@ def read_ms_i(name, data_column, field, spectral_window, pol_indices, pol_summat # Check if data column is available and get shape with ms_table(name) as t: nmspol = t.getcol(data_column, startrow=0, nrow=3).shape[2] - nrow = t.nrow() + nrow = t.nrows() print("Measurement set visibilities:") print(f" shape: ({nrow}, {_ms_nchannels(name, spectral_window)}, {nmspol})") @@ -180,10 +180,10 @@ def read_ms_i(name, data_column, field, spectral_window, pol_indices, pol_summat if pol_summation: npol = 1 else: - if pol_indices != slice(None): - npol = len(pol_indices) + if pol_indices == slice(None): + npol = nmspol else: - npol = npol + npol = len(pol_indices) shp = (nrealrows, nrealchan, npol) vis = np.empty(shp, dtype=np.complex64) wgt = np.empty(shp, dtype=np.float32) -- GitLab From c68d914cd750f0fad58d97630e35a5263acc74f0 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Tue, 31 Aug 2021 13:32:35 +0200 Subject: [PATCH 27/89] Cosmetics --- resolve/data/ms_import.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/resolve/data/ms_import.py b/resolve/data/ms_import.py index 0678a31c..f16ce733 100644 --- a/resolve/data/ms_import.py +++ b/resolve/data/ms_import.py @@ -319,9 +319,7 @@ def _first_pass(ms, field, spectral_window, channels, pol_indices, pol_summation assert tflags.ndim == 3 tflags = np.all(tflags, axis=2) active_rows[start:stop] = np.invert(np.all(tflags, axis=1)) - active_channels = np.logical_or( - active_channels, np.invert(np.all(tflags, axis=0)) - ) + active_channels = np.logical_or(active_channels, np.invert(np.all(tflags, axis=0))) start = stop return active_rows, active_channels -- GitLab From 033b793cfff097dbe6b79d61c9166f5a64bb66bd Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Tue, 31 Aug 2021 17:16:52 +0200 Subject: [PATCH 28/89] Fixup --- resolve/data/ms_import.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/resolve/data/ms_import.py b/resolve/data/ms_import.py index f16ce733..b4abbefb 100644 --- a/resolve/data/ms_import.py +++ b/resolve/data/ms_import.py @@ -143,6 +143,13 @@ def _determine_weighting(ms): if "WEIGHT_SPECTRUM" in t.colnames(): weightcol = "WEIGHT_SPECTRUM" fullwgt = True + + # TODO Find a better way to detect whether weight spectrum is available + try: + t.getcol("WEIGHT_SPECTRUM", nrow=1) + except RuntimeError: + weightcol = "WEIGHT" + fullwgt = False else: weightcol = "WEIGHT" fullwgt = False -- GitLab From 16fe53a4b62015e2832be7cd095cb5f53afd20d7 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Tue, 31 Aug 2021 17:49:46 +0200 Subject: [PATCH 29/89] Alma beam: return beam if cache is empty --- resolve/library/primary_beams.py | 1 + 1 file changed, 1 insertion(+) diff --git a/resolve/library/primary_beams.py b/resolve/library/primary_beams.py index e43abe66..f351d55a 100644 --- a/resolve/library/primary_beams.py +++ b/resolve/library/primary_beams.py @@ -241,6 +241,7 @@ def alma_beam_func(D, d, freq, x, use_cache=False): except FileNotFoundError: arr = _compute_alma_beam(D, d, freq, x) np.save(fname, arr) + return arr def _compute_alma_beam(D, d, freq, x): -- GitLab From 1ac72769462dcae96933ef5a5e9cdeaf3df2fd1d Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Tue, 31 Aug 2021 18:16:18 +0200 Subject: [PATCH 30/89] Warning concerning POINTING table --- resolve/data/ms_import.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/resolve/data/ms_import.py b/resolve/data/ms_import.py index b4abbefb..6d119e6b 100644 --- a/resolve/data/ms_import.py +++ b/resolve/data/ms_import.py @@ -124,8 +124,7 @@ def ms2observations(ms, data_column, with_calib_info, spectral_window, print(f"{ms}, field #{ifield} is empty or fully flagged") observations.append(None) continue - if mm["ptg"] is not None: - raise NotImplementedError + print("WARNING: POINTING table is ignored for now") antpos = AntennaPositions(mm["uvw"], mm["ant1"], mm["ant2"], mm["time"]) obs = Observation(antpos, mm["vis"], mm["wgt"], polobj, mm["freq"], auxiliary_tables=auxtables) -- GitLab From 660fc09fbdfca0ad61005cedd7b5029067ebd239 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Tue, 31 Aug 2021 18:36:37 +0200 Subject: [PATCH 31/89] Add source name to Observation --- resolve/data/observation.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/resolve/data/observation.py b/resolve/data/observation.py index 73407d08..f5efa47a 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -497,6 +497,12 @@ class Observation(BaseObservation): def auxiliary_table(self, name): return self._auxiliary_tables[name] + @property + def source_name(self): + if "FIELD" in self._auxiliary_tables: + return self._auxiliary_tables["FIELD"]["NAME"][0] + raise NotImplementedError("FIELD subtable not available.") + def effective_uvw(self): out = np.einsum("ij,k->jik", self.uvw, self._freq / SPEEDOFLIGHT) my_asserteq(out.shape, (3, self.nrow, self.nfreq)) -- GitLab From 0e67047a7ecd39b5214269ddeb731b39f343009b Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Tue, 31 Aug 2021 19:48:24 +0200 Subject: [PATCH 32/89] plotter.py: add basic fits support --- resolve/plotter.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/resolve/plotter.py b/resolve/plotter.py index 9d4de344..90c09dcc 100644 --- a/resolve/plotter.py +++ b/resolve/plotter.py @@ -27,10 +27,16 @@ def makedirs(*args, **kwargs): class Plotter: def __init__(self, fileformat, directory): self._nifty, self._calib, self._hist = [], [], [] + self._fits = [] self._f = fileformat self._dir = directory makedirs(self._dir, exist_ok=True) + @onlymaster + def add_fits(self, name, operator, direction=None): + self._fits.append({"operator": operator, "title": str(name), + "direction": direction}) + @onlymaster def add_npy(self, name, operator, nsamples=0): """Save output of operator as npy files. @@ -92,6 +98,13 @@ class Plotter: if fname is None: continue _plot_histograms(state, fname, 20, 0.5, postop=op) + for ii, obj in enumerate(self._fits): + from .fits import field2fits + op, fname = self._plot_init(obj, state, identifier) + sc = ift.StatCalculator() + for ss in state: + sc.add(op.force(ss)) + field2fits(sc.mean, f"{fname}.fits", True, obj["direction"]) mydir = join(self._dir, "zzz_latent") makedirs(mydir, exist_ok=True) _plot_histograms(state, join(mydir, f"{identifier}.{self._f}"), 5, 0.5) -- GitLab From f55f8c5d33f03b1c356e4f8a06bd99dc159a218f Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Wed, 1 Sep 2021 13:51:03 +0200 Subject: [PATCH 33/89] Add polarization space --- resolve/__init__.py | 1 + resolve/polarization_space.py | 45 +++++++++++++++++++++++++++++++++++ 2 files changed, 46 insertions(+) create mode 100644 resolve/polarization_space.py diff --git a/resolve/__init__.py b/resolve/__init__.py index 4a1ade5b..68cdefef 100644 --- a/resolve/__init__.py +++ b/resolve/__init__.py @@ -23,6 +23,7 @@ from .plotter import MfPlotter, Plotter from .points import PointInserter from .data.polarization import Polarization from .polarization_matrix_exponential import * +from .polarization_space import * from .response import MfResponse, ResponseDistributor, StokesIResponse, SingleResponse from .simple_operators import * from .util import * diff --git a/resolve/polarization_space.py b/resolve/polarization_space.py new file mode 100644 index 00000000..380681fb --- /dev/null +++ b/resolve/polarization_space.py @@ -0,0 +1,45 @@ +# This program 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. +# +# This program 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 this program. If not, see <http://www.gnu.org/licenses/>. +# +# Copyright(C) 2013-2021 Max-Planck-Society + +import nifty8 as ift + +from .util import my_assert + + +class PolarizationSpace(ift.UnstructuredDomain): + """ + + Parameters + ---------- + coordinates : np.ndarray + Must be sorted and strictly ascending. + """ + + _needed_for_hash = ["_shape", "_lbl"] + + def __init__(self, polarization_labels): + if isinstance(polarization_labels, str): + polarization_labels = (polarization_labels,) + self._lbl = tuple(polarization_labels) + for ll in self._lbl: + my_assert(ll in ["I", "Q", "U", "V", "LL", "LR", "RL", "RR", "XX", "XY", "YX", "YY"]) + super(PolarizationSpace, self).__init__(len(self._lbl)) + + def __repr__(self): + return f"PolarizationSpace({self._lbl})" + + @property + def labels(self): + return self._lbl -- GitLab From e3577ea3ae45f6152e83ec442999f7ca669d8dfc Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Wed, 1 Sep 2021 13:53:39 +0200 Subject: [PATCH 34/89] Observation: add ant and time to interface --- resolve/data/observation.py | 12 ++++++++++++ resolve/polarization_space.py | 17 +++-------------- 2 files changed, 15 insertions(+), 14 deletions(-) diff --git a/resolve/data/observation.py b/resolve/data/observation.py index 73407d08..269d00f8 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -485,6 +485,18 @@ class Observation(BaseObservation): def antenna_positions(self): return self._antpos + @property + def ant1(self): + return self._antpos.ant1 + + @property + def ant2(self): + return self._antpos.ant2 + + @property + def time(self): + return self._antpos.time + @property def direction(self): if self._auxiliary_tables is not None and "FIELD" in self._auxiliary_tables: diff --git a/resolve/polarization_space.py b/resolve/polarization_space.py index 380681fb..5cb209e8 100644 --- a/resolve/polarization_space.py +++ b/resolve/polarization_space.py @@ -1,17 +1,6 @@ -# This program 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. -# -# This program 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 this program. If not, see <http://www.gnu.org/licenses/>. -# -# Copyright(C) 2013-2021 Max-Planck-Society +# SPDX-License-Identifier: GPL-3.0-or-later +# Copyright(C) 2021 Max-Planck-Society +# Author: Philipp Arras, Jakob Knollmüller import nifty8 as ift -- GitLab From 6a0ead301e59f93e7b027848ee91fc756c7d6ba3 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Wed, 1 Sep 2021 15:58:18 +0200 Subject: [PATCH 35/89] SingleResponse: make mask optional --- resolve/data/observation.py | 8 ++++++++ resolve/response.py | 5 +++-- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/resolve/data/observation.py b/resolve/data/observation.py index 269d00f8..ba72f0ad 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -452,6 +452,14 @@ class Observation(BaseObservation): self._antpos, vis, wgt, Polarization.trivial(), self._freq, self._auxiliary_tables ) + def restrict_by_time(self, tmin, tmax): + ind = np.logical_and(tmin <= self.time , self.time < tmax) + return self[ind] + + def restrict_by_freq(self, fmin, fmax): + ind = np.logical_and(fmin <= self.freq, self.freq < fmax) + return self.get_freqs_by_slice(ind) + def restrict_to_stokesi(self): # FIXME Do I need to change something in self._auxiliary_tables? my_asserteq(self._vis.shape[0], 4) diff --git a/resolve/response.py b/resolve/response.py index c9d66ee3..5e074cb9 100644 --- a/resolve/response.py +++ b/resolve/response.py @@ -222,7 +222,7 @@ class FullResponse(ift.LinearOperator): class SingleResponse(ift.LinearOperator): - def __init__(self, domain, uvw, freq, mask, facets=(1, 1)): + def __init__(self, domain, uvw, freq, mask=None, facets=(1, 1)): my_assert_isinstance(facets, tuple) for ii in range(2): if domain.shape[ii] % facets[ii] != 0: @@ -234,7 +234,6 @@ class SingleResponse(ift.LinearOperator): self._args = { "uvw": uvw, "freq": freq, - "mask": mask.astype(np.uint8), "pixsize_x": self._domain[0].distances[0], "pixsize_y": self._domain[0].distances[1], "epsilon": epsilon(), @@ -242,6 +241,8 @@ class SingleResponse(ift.LinearOperator): "nthreads": nthreads(), "flip_v": True, } + if mask is not None: + self._args["mask"] = mask.astype(np.uint8) self._vol = self._domain[0].scalar_dvol self._target_dtype = np_dtype(True) self._domain_dtype = np_dtype(False) -- GitLab From 9096a07d9254ec2bf29a992ee4ce03a7e87e5b84 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Wed, 1 Sep 2021 15:58:43 +0200 Subject: [PATCH 36/89] Observation: restrict by time/freq via searchsorted --- resolve/data/observation.py | 26 +++++++++++++++++++------- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/resolve/data/observation.py b/resolve/data/observation.py index ba72f0ad..b194b235 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -294,6 +294,8 @@ class Observation(BaseObservation): my_assert(np.all(np.isfinite(vis[weight > 0.]))) my_assert(np.all(np.isfinite(weight))) + my_assert(np.all(np.diff(freq))) + vis.flags.writeable = False weight.flags.writeable = False @@ -452,13 +454,23 @@ class Observation(BaseObservation): self._antpos, vis, wgt, Polarization.trivial(), self._freq, self._auxiliary_tables ) - def restrict_by_time(self, tmin, tmax): - ind = np.logical_and(tmin <= self.time , self.time < tmax) - return self[ind] - - def restrict_by_freq(self, fmin, fmax): - ind = np.logical_and(fmin <= self.freq, self.freq < fmax) - return self.get_freqs_by_slice(ind) + def restrict_by_time(self, tmin, tmax, with_index=False): + my_assert(all(np.diff(self.time) >= 0)) + start, stop = np.searchsorted(self.time, [tmin, tmax]) + ind = slice(start, stop) + res = self[ind] + if with_index: + return res, ind + return res + + def restrict_by_freq(self, fmin, fmax, with_index=False): + my_assert(all(np.diff(self.freq) > 0)) + start, stop = np.searchsorted(self.freq, [fmin, fmax]) + ind = slice(start, stop) + res = self.get_freqs_by_slice(ind) + if with_index: + return res, ind + return res def restrict_to_stokesi(self): # FIXME Do I need to change something in self._auxiliary_tables? -- GitLab From d3e529770efa462c68ca0c50b384bbdd3bfb42e2 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Wed, 1 Sep 2021 15:59:12 +0200 Subject: [PATCH 37/89] Cosmetics --- resolve/polarization_space.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/resolve/polarization_space.py b/resolve/polarization_space.py index 5cb209e8..dad3e786 100644 --- a/resolve/polarization_space.py +++ b/resolve/polarization_space.py @@ -27,7 +27,7 @@ class PolarizationSpace(ift.UnstructuredDomain): super(PolarizationSpace, self).__init__(len(self._lbl)) def __repr__(self): - return f"PolarizationSpace({self._lbl})" + return f"PolarizationSpace(polarization_labels={self._lbl})" @property def labels(self): -- GitLab From 13be61630acf96b494562ee65b644c0bc7637b0b Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 2 Sep 2021 10:40:12 +0200 Subject: [PATCH 38/89] Cleanup --- resolve/multi_frequency/irg_space.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/resolve/multi_frequency/irg_space.py b/resolve/multi_frequency/irg_space.py index 70a07f71..a2a23bc7 100644 --- a/resolve/multi_frequency/irg_space.py +++ b/resolve/multi_frequency/irg_space.py @@ -12,8 +12,6 @@ # along with this program. If not, see <http://www.gnu.org/licenses/>. # # Copyright(C) 2013-2020 Max-Planck-Society -# -# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. import numpy as np -- GitLab From deb6e5bf5b081a906a15d5588721a3c4acb0f270 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 3 Sep 2021 18:06:43 +0200 Subject: [PATCH 39/89] Observation: Use PolarizationSpace for data domain --- resolve/data/observation.py | 3 ++- resolve/data/polarization.py | 9 +++++++++ resolve/likelihood.py | 2 +- resolve/multi_frequency/operators.py | 6 +++--- resolve/response.py | 8 +++++--- test/test_general.py | 12 +++++------- 6 files changed, 25 insertions(+), 15 deletions(-) diff --git a/resolve/data/observation.py b/resolve/data/observation.py index b194b235..6fb42d80 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -18,7 +18,8 @@ from .polarization import Polarization class BaseObservation: @property def _dom(self): - dom = [ift.UnstructuredDomain(ss) for ss in self._vis.shape] + pol_dom = self.polarization.space + dom = [pol_dom] + [ift.UnstructuredDomain(ss) for ss in self._vis.shape[1:]] return ift.makeDomain(dom) @property diff --git a/resolve/data/polarization.py b/resolve/data/polarization.py index 3bbbd199..48b2c79c 100644 --- a/resolve/data/polarization.py +++ b/resolve/data/polarization.py @@ -68,6 +68,15 @@ class Polarization: def to_str_list(self): return [TABLE[ii] for ii in self._ind] + @property + def space(self): + from ..polarization_space import PolarizationSpace + if self == Polarization.trivial(): + x = "I" + else: + x = [TABLE[ii] for ii in self._ind] + return PolarizationSpace(x) + @staticmethod def from_list(lst): return Polarization(lst) diff --git a/resolve/likelihood.py b/resolve/likelihood.py index 7e75481a..54c9ed46 100644 --- a/resolve/likelihood.py +++ b/resolve/likelihood.py @@ -1,5 +1,5 @@ # SPDX-License-Identifier: GPL-3.0-or-later -# Copyright(C) 2020 Max-Planck-Society +# Copyright(C) 2020-2021 Max-Planck-Society # Author: Philipp Arras from functools import reduce diff --git a/resolve/multi_frequency/operators.py b/resolve/multi_frequency/operators.py index 11aae3db..b77f98b2 100644 --- a/resolve/multi_frequency/operators.py +++ b/resolve/multi_frequency/operators.py @@ -2,10 +2,10 @@ # Copyright(C) 2013-2021 Max-Planck-Society # Authors: Philipp Frank, Philipp Arras, Philipp Haim -import numpy as np - import nifty8 as ift +import numpy as np +from ..polarization_space import PolarizationSpace from ..util import my_assert, my_assert_isinstance, my_asserteq from .irg_space import IRGSpace @@ -91,7 +91,7 @@ class MfWeightingInterpolation(ift.LinearOperator): my_asserteq(domain.shape[0], eff_uvw.shape[2]) # freqaxis self._domain = domain nrow, nfreq = eff_uvw.shape[1:] - tgt = [ift.UnstructuredDomain(aa) for aa in [1, nrow, nfreq]] + tgt = [PolarizationSpace("I")] + [ift.UnstructuredDomain(aa) for aa in [nrow, nfreq]] self._target = ift.DomainTuple.make(tgt) self._capability = self.TIMES | self.ADJOINT_TIMES # FIXME Try to unify all those operators which loop over freq dimension diff --git a/resolve/response.py b/resolve/response.py index 5e074cb9..fcb097b3 100644 --- a/resolve/response.py +++ b/resolve/response.py @@ -41,7 +41,7 @@ def StokesIResponse(observation, domain): return contr.adjoint @ sr0 elif npol == 2: sr1 = SingleResponse(domain, observation.uvw, observation.freq, mask[1]) - return ResponseDistributor(sr0, sr1) + return ResponseDistributor(observation.polarization.space, sr0, sr1) raise RuntimeError @@ -186,7 +186,7 @@ class MfResponse(ift.LinearOperator): class ResponseDistributor(ift.LinearOperator): - def __init__(self, *ops): + def __init__(self, iterating_target_space, *ops): dom, tgt = ops[0].domain, ops[0].target cap = self.TIMES | self.ADJOINT_TIMES for op in ops: @@ -194,8 +194,10 @@ class ResponseDistributor(ift.LinearOperator): my_assert(dom is op.domain) my_assert(tgt is op.target) my_assert(self.TIMES & op.capability, self.ADJOINT_TIMES & op.capability) + my_asserteq(len(iterating_target_space.shape), 1) + my_asserteq(len(ops), iterating_target_space.shape[0]) self._domain = ift.makeDomain(dom) - self._target = ift.makeDomain((ift.UnstructuredDomain(len(ops)), *tgt)) + self._target = ift.makeDomain((iterating_target_space, *tgt)) self._capability = cap self._ops = ops diff --git a/test/test_general.py b/test/test_general.py index 13845c50..4c91daf9 100644 --- a/test/test_general.py +++ b/test/test_general.py @@ -1,5 +1,5 @@ # SPDX-License-Identifier: GPL-3.0-or-later -# Copyright(C) 2020 Max-Planck-Society +# Copyright(C) 2020-2021 Max-Planck-Society # Author: Philipp Arras from os.path import join @@ -146,7 +146,7 @@ def test_calibration_likelihood(time_mode): npol, nfreq = obs[0].npol, obs[0].nfreq total_N = npol * nants * nfreq dom = [ - ift.UnstructuredDomain(npol), + obs[0].polarization.space, ift.UnstructuredDomain(len(uants)), time_domain, ift.UnstructuredDomain(nfreq), @@ -213,10 +213,8 @@ def test_calibration_distributor(obs): tgt = obs.vis.domain utimes = rve.unique_times(obs) uants = obs.antenna_positions.unique_antennas() - dom = [ - ift.UnstructuredDomain(nn) - for nn in [obs.npol, len(uants), len(utimes), obs.nfreq] - ] + dom = [obs.polarization.space] + \ + [ift.UnstructuredDomain(nn) for nn in [len(uants), len(utimes), obs.nfreq]] uants = rve.unique_antennas(obs) time_dct = {aa: ii for ii, aa in enumerate(utimes)} antenna_dct = {aa: ii for ii, aa in enumerate(uants)} @@ -259,7 +257,7 @@ def test_response_distributor(): dom = ift.UnstructuredDomain(2), ift.UnstructuredDomain(4) op0 = ift.makeOp(ift.makeField(dom, np.arange(8).reshape(2, -1))) op1 = ift.makeOp(ift.makeField(dom, 2 * np.arange(8).reshape(2, -1))) - op = rve.response.ResponseDistributor(op0, op1) + op = rve.response.ResponseDistributor(ift.UnstructuredDomain(2), op0, op1) ift.extra.check_linear_operator(op) -- GitLab From f2c9eb723eb9b92be8700760bc1daf7a8c2409fd Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Mon, 6 Sep 2021 12:13:13 +0200 Subject: [PATCH 40/89] Data import: Fixup polarization selection --- resolve/data/ms_import.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/resolve/data/ms_import.py b/resolve/data/ms_import.py index 6d119e6b..45c7e837 100644 --- a/resolve/data/ms_import.py +++ b/resolve/data/ms_import.py @@ -100,8 +100,8 @@ def ms2observations(ms, data_column, with_calib_info, spectral_window, else: if polarizations == "stokesi": polarizations = ["LL", "RR"] if polobj.circular() else ["XX", "YY"] - polobj = polobj.restrict_by_name(polarizations) pol_ind = [polobj.to_str_list().index(ii) for ii in polarizations] + polobj = polobj.restrict_by_name(polarizations) pol_summation = False observations = [] -- GitLab From 0aec963e012bcbf9a8a6bcb2a229c1fd8a9cb467 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Tue, 7 Sep 2021 13:28:41 +0200 Subject: [PATCH 41/89] Add station names --- resolve/data/observation.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/resolve/data/observation.py b/resolve/data/observation.py index f5efa47a..5fef0166 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -503,6 +503,15 @@ class Observation(BaseObservation): return self._auxiliary_tables["FIELD"]["NAME"][0] raise NotImplementedError("FIELD subtable not available.") + @property + def station_names(self): + """The index of the resulting list is the same as in self.ant1 or self.ant2.""" + if "ANTENNA" in self._auxiliary_tables: + tab = self._auxiliary_tables["ANTENNA"] + return [f"{a} {b}" for a, b in zip(tab["NAME"], tab["STATION"])] + raise NotImplementedError("ANTENNA subtable not available.") + + def effective_uvw(self): out = np.einsum("ij,k->jik", self.uvw, self._freq / SPEEDOFLIGHT) my_asserteq(out.shape, (3, self.nrow, self.nfreq)) -- GitLab From 6998661ff25dc7093c1790316c19253594358146 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Sun, 12 Sep 2021 11:39:33 +0200 Subject: [PATCH 42/89] Minor --- resolve/calibration.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/resolve/calibration.py b/resolve/calibration.py index 0fb9b8fb..bcd322c8 100644 --- a/resolve/calibration.py +++ b/resolve/calibration.py @@ -37,6 +37,8 @@ class CalibrationDistributor(ift.LinearOperator): my_assert(np.max(ant_col) < len(antenna_dct)) self._domain = ift.DomainTuple.make(domain) self._target = ift.DomainTuple.make(target) + my_asserteq(len(self._domain), 4) + my_asserteq(len(self._target), 3) self._capability = self.TIMES | self.ADJOINT_TIMES self._nantennas = len(antenna_dct) @@ -101,9 +103,7 @@ class MyLinearInterpolator(ift.LinearOperator): my_assert(np.min(time_col) >= 0) my_assert(np.max(time_col) < self._domain.shape[1] * dom[0].distances[1]) my_assert(np.issubdtype(ant_col.dtype, np.integer)) - my_assert( - np.issubdtype(time_col.dtype, np.floating if floattime else np.integer) - ) + my_assert(np.issubdtype(time_col.dtype, np.floating if floattime else np.integer)) def apply(self, x, mode): self._check_input(x, mode) -- GitLab From 52dbf9a068803f02b3a6fc69f381b0231b088abb Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 16 Sep 2021 15:27:45 +0200 Subject: [PATCH 43/89] Observation: Add legacy_load --- resolve/data/observation.py | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/resolve/data/observation.py b/resolve/data/observation.py index 5fef0166..fd44cde7 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -369,6 +369,30 @@ class Observation(BaseObservation): auxtables ) + @staticmethod + def legacy_load(file_name, lo_hi_index=None): + """Provide potentially incomplete interface for loading legacy npz files.""" + dct = np.load(file_name) + antpos = [] + for ii in range(4): + val = dct[f"antpos{ii}"] + if val.size == 0: + val = None + antpos.append(val) + antpos = AntennaPositions.from_list(antpos) + pol = Polarization.from_list(dct["polarization"]) + vis = dct["vis"] + wgt = dct["weight"] + freq = dct["freq"] + if lo_hi_index is not None: + slc = slice(*lo_hi_index) + # Convert view into its own array + vis = vis[..., slc].copy() + wgt = wgt[..., slc].copy() + freq = freq[slc].copy() + del dct + return Observation(antpos, vis, wgt, pol, freq) + def flags_to_nan(self): if self.fraction_useful == 1.: return self @@ -511,7 +535,6 @@ class Observation(BaseObservation): return [f"{a} {b}" for a, b in zip(tab["NAME"], tab["STATION"])] raise NotImplementedError("ANTENNA subtable not available.") - def effective_uvw(self): out = np.einsum("ij,k->jik", self.uvw, self._freq / SPEEDOFLIGHT) my_asserteq(out.shape, (3, self.nrow, self.nfreq)) -- GitLab From 3d78177fdd0b4a0832325bc07b15ef43e3db6631 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 16 Sep 2021 15:34:47 +0200 Subject: [PATCH 44/89] Observation: add legacy load --- test/test_general.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/test_general.py b/test/test_general.py index 13845c50..44e255b9 100644 --- a/test/test_general.py +++ b/test/test_general.py @@ -67,6 +67,10 @@ def test_sliced_readin(slc): np.testing.assert_equal(ch0, ch) +def test_legacy_load(): + rve.Observation.legacy_load(f"{direc}legacy.npz") + + def try_operator(op): pos = ift.from_random(op.domain) op(pos) -- GitLab From 762b7da1931d3099abd77c97a25ad13070c2d462 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 17 Sep 2021 15:14:27 +0200 Subject: [PATCH 45/89] Observation: add time_average --- resolve/data/observation.py | 47 +++++++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/resolve/data/observation.py b/resolve/data/observation.py index 6fb42d80..21d3bc48 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -498,6 +498,53 @@ class Observation(BaseObservation): self._auxiliary_tables ) + def time_average(self, list_of_timebins): + # FIXME check that timebins do not overlap + # time, ant1, ant2 + ts = self._antpos.time + row_to_bin_map = np.empty(ts.shape) + row_to_bin_map[:] = np.nan + + for ii, (lo, hi) in enumerate(list_of_timebins): + ind = np.logical_and(ts >= lo, ts < hi) + assert np.all(np.isnan(row_to_bin_map[ind])) + row_to_bin_map[ind] = ii + + assert np.all(~np.isnan(row_to_bin_map)) + row_to_bin_map = row_to_bin_map.astype(int) + + ant1 = self._antpos.ant1 + ant2 = self._antpos.ant2 + assert np.max(ant1) < 1000 + assert np.max(ant2) < 1000 + assert np.max(row_to_bin_map) < np.iinfo(np.dtype("int64")).max / 1000000 + atset = set(zip(ant1, ant2, row_to_bin_map)) + dct = {aa: ii for ii, aa in enumerate(atset)} + dct_inv = {yy: xx for xx, yy in dct.items()} + masterindex = np.array([dct[(a1, a2, tt)] for a1, a2, tt in zip(ant1, ant2, row_to_bin_map)]) + + vis, wgt = self.vis.val, self.weight.val + new_vis = np.empty((self.npol, len(atset), self.nfreq), dtype=self.vis.dtype) + new_wgt = np.empty((self.npol, len(atset), self.nfreq), dtype=self.weight.dtype) + for pol in range(self.npol): + for freq in range(self.nfreq): + enum = np.bincount(masterindex, weights=vis[pol, :, freq].real*wgt[pol, :, freq]) + enum = enum + 1j*np.bincount(masterindex, weights=vis[pol, :, freq].imag*wgt[pol, :, freq]) + denom = np.bincount(masterindex, weights=wgt[pol, :, freq]) + new_vis[pol, :, freq] = enum/denom + new_wgt[pol, :, freq] = denom + + new_times = np.array([dct_inv[ii][2] for ii in range(len(atset))]) + new_times = np.mean(np.array(list_of_timebins), axis=1)[new_times] + + new_ant1 = np.array([dct_inv[ii][0] for ii in range(len(atset))]) + new_ant2 = np.array([dct_inv[ii][1] for ii in range(len(atset))]) + new_uvw = ap.uvw + # FIXME Find correct uvw + ap = self._antpos + ap = AntennaPositions(new_uvw, new_ant1, new_ant2, new_times) + return Observation(ap, new_vis, new_wgt, self._polarization, self._freq, self._auxiliary_tables) + @property def uvw(self): return self._antpos.uvw -- GitLab From 89689d1dc8dd7c126609e110aa259f7a0dbbabe6 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Wed, 27 Oct 2021 10:09:29 +0200 Subject: [PATCH 46/89] Observation: average of uvw for time binning --- resolve/data/observation.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/resolve/data/observation.py b/resolve/data/observation.py index e5636b54..a290b499 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -556,15 +556,19 @@ class Observation(BaseObservation): new_vis[pol, :, freq] = enum/denom new_wgt[pol, :, freq] = denom + new_uvw = np.empty((len(atset), 3), dtype=self._antpos.uvw.dtype) + new_uvw[()] = np.nan + denom = np.bincount(masterindex) + # Assumption: Uvw value for averaged data is average of uvw values of finely binned data + for ii in range(3): + new_uvw[:, ii] = np.bincount(masterindex, weights=self._antpos.uvw[:, ii]) / denom + assert np.sum(np.isnan(new_uvw)) == 0 + new_times = np.array([dct_inv[ii][2] for ii in range(len(atset))]) new_times = np.mean(np.array(list_of_timebins), axis=1)[new_times] new_ant1 = np.array([dct_inv[ii][0] for ii in range(len(atset))]) new_ant2 = np.array([dct_inv[ii][1] for ii in range(len(atset))]) - new_uvw = ap.uvw - raise NotImplementedError("FIXME The implementation of this is not complete.") - # FIXME Find correct uvw - ap = self._antpos ap = AntennaPositions(new_uvw, new_ant1, new_ant2, new_times) return Observation(ap, new_vis, new_wgt, self._polarization, self._freq, self._auxiliary_tables) -- GitLab From dd48000dcceddabae025489cb5b901b454a6605d Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Wed, 27 Oct 2021 10:40:08 +0200 Subject: [PATCH 47/89] Adjust to new nifty interface --- resolve/likelihood.py | 2 +- resolve/minimization.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/resolve/likelihood.py b/resolve/likelihood.py index 54c9ed46..642c7159 100644 --- a/resolve/likelihood.py +++ b/resolve/likelihood.py @@ -55,7 +55,7 @@ def _build_gauss_lh_nres(op, mean, invcov): my_assert_isinstance(mean, invcov, (ift.Field, ift.MultiField)) my_asserteq(op.target, mean.domain, invcov.domain) - lh = ift.GaussianEnergy(mean=mean, inverse_covariance=ift.makeOp(invcov)) @ op + lh = ift.GaussianEnergy(mean=mean, inverse_covariance=ift.makeOp(invcov, sampling_dtype=mean.dtype)) @ op return _Likelihood(lh, mean, lambda x: ift.makeOp(invcov), op) diff --git a/resolve/minimization.py b/resolve/minimization.py index 32ef0a72..dbf004c5 100644 --- a/resolve/minimization.py +++ b/resolve/minimization.py @@ -28,7 +28,7 @@ class Minimization: else: my_assert(n_samples > 0) dct = { - "mean": position, + "position": position, "hamiltonian": operator, "n_samples": n_samples, "minimizer_sampling": None, @@ -38,7 +38,7 @@ class Minimization: "comm": comm, "nanisinf": True, } - self._e = ift.SampledKL(**dct) + self._e = ift.SampledKLEnergy(**dct) self._n, self._m = dct["n_samples"], dct["mirror_samples"] def minimize(self, minimizer): -- GitLab From ab7748eec21e57b4723c8e3b478e20cb14a0ceb2 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 28 Oct 2021 00:16:36 +0200 Subject: [PATCH 48/89] Observation: properly sort time when binning --- resolve/data/observation.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/resolve/data/observation.py b/resolve/data/observation.py index a290b499..bc202224 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -531,6 +531,7 @@ class Observation(BaseObservation): ind = np.logical_and(ts >= lo, ts < hi) assert np.all(np.isnan(row_to_bin_map[ind])) row_to_bin_map[ind] = ii + assert np.all(np.diff(row_to_bin_map) >= 0) assert np.all(~np.isnan(row_to_bin_map)) row_to_bin_map = row_to_bin_map.astype(int) @@ -540,7 +541,9 @@ class Observation(BaseObservation): assert np.max(ant1) < 1000 assert np.max(ant2) < 1000 assert np.max(row_to_bin_map) < np.iinfo(np.dtype("int64")).max / 1000000 - atset = set(zip(ant1, ant2, row_to_bin_map)) + atset = np.array(list(set(zip(ant1, ant2, row_to_bin_map)))) + atset = atset[np.lexsort(atset.T)] + atset = tuple(map(tuple, atset)) dct = {aa: ii for ii, aa in enumerate(atset)} dct_inv = {yy: xx for xx, yy in dct.items()} masterindex = np.array([dct[(a1, a2, tt)] for a1, a2, tt in zip(ant1, ant2, row_to_bin_map)]) @@ -567,6 +570,8 @@ class Observation(BaseObservation): new_times = np.array([dct_inv[ii][2] for ii in range(len(atset))]) new_times = np.mean(np.array(list_of_timebins), axis=1)[new_times] + assert np.all(np.diff(new_times) >= 0) + new_ant1 = np.array([dct_inv[ii][0] for ii in range(len(atset))]) new_ant2 = np.array([dct_inv[ii][1] for ii in range(len(atset))]) ap = AntennaPositions(new_uvw, new_ant1, new_ant2, new_times) -- GitLab From 13f164ec7b06a24ff357d35e3c1fca8e12e8232e Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 29 Oct 2021 11:43:28 +0200 Subject: [PATCH 49/89] Add VLBI averaging --- resolve/__init__.py | 1 + resolve/data/averaging.py | 68 +++++++++++++++++++++++++++++++++++++ resolve/data/observation.py | 9 ++--- 3 files changed, 74 insertions(+), 4 deletions(-) create mode 100644 resolve/data/averaging.py diff --git a/resolve/__init__.py b/resolve/__init__.py index 49833d48..58bb5754 100644 --- a/resolve/__init__.py +++ b/resolve/__init__.py @@ -4,6 +4,7 @@ from .library.calibrators import * from .library.primary_beams import * from .constants import * from .data.direction import * +from .data.averaging import * from .fits import field2fits, fits2field from .global_config import * from .likelihood import * diff --git a/resolve/data/averaging.py b/resolve/data/averaging.py new file mode 100644 index 00000000..65eb7a5e --- /dev/null +++ b/resolve/data/averaging.py @@ -0,0 +1,68 @@ +# SPDX-License-Identifier: GPL-3.0-or-later +# Copyright(C) 2019-2021 Max-Planck-Society +# Author: Martin Reinecke + +import numpy as np + + +def fair_share_averaging(ts_per_bin, times, gap_time): + """Bin time stamps fairly. + + Parameters + ---------- + ts_per_bin : int + Minimum number of data points in one average. + times : np.array + One-dimensional array of time stamps. + gap_time : float + Maximum time difference that can occur in one average. + """ + if ts_per_bin == 1: + return times + if ts_per_bin < 1: + raise ValueError + + times = np.sort(times.copy()) + tbins = [] + nval = len(times) + i0 = 0 + + # The following algorithm was designed to write the same time stamp for a + # given bin into the original array. Now we need only the bounds of the time + # bins. Thus, the following implementation could be improved to directly + # compute these. + while i0 < nval: + i = i0+1 + nscan = 1 # number of different time stamps in the scan + while i < nval and times[i]-times[i-1] < gap_time: # as long as there are less than x seconds between time stamps, we assume we are in the same scan + if times[i] != times[i-1]: + nscan += 1 + i += 1 + nbin = max(1, nscan//ts_per_bin) # how many bins to use for this scan + for j in range(nbin): + n = _fair_share(nscan, nbin, j) + i = i0+1 + icnt = 0 + oldtime = times[i0] + while i < nval and icnt < n: + if times[i] != oldtime: + icnt += 1 + oldtime = times[i] + if icnt < n: + if icnt == n-1: + tbins += [(times[i0], times[i], times[i]-times[i0])] + times[i] = times[i0] # give all values in this bin the time stamp of the first value + i += 1 + i0 = i + tbsize = np.array([t[2] for t in tbins]) + print("Size time bins:") + print(f" min: {np.min(tbsize):.1f}s\n max: {np.max(tbsize):.1f}s") + print(f" mean: {np.mean(tbsize):.1f}s\n median: {np.median(tbsize):.1f}s") + times = np.unique(times) + times = np.hstack([times, np.inf]) + time_bins = np.array([times[:-1], times[1:]]).T + return time_bins + + +def _fair_share(n, nshare, ishare): + return n//nshare + (ishare < (n % nshare)) diff --git a/resolve/data/observation.py b/resolve/data/observation.py index bc202224..7c41e868 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -560,16 +560,17 @@ class Observation(BaseObservation): new_wgt[pol, :, freq] = denom new_uvw = np.empty((len(atset), 3), dtype=self._antpos.uvw.dtype) - new_uvw[()] = np.nan + new_times = np.empty(len(atset), dtype=self._antpos.time.dtype) + new_uvw[()] = new_times[()] = np.nan denom = np.bincount(masterindex) # Assumption: Uvw value for averaged data is average of uvw values of finely binned data for ii in range(3): new_uvw[:, ii] = np.bincount(masterindex, weights=self._antpos.uvw[:, ii]) / denom + new_times = np.bincount(row_to_bin_map, weights=self._antpos.time) / np.bincount(row_to_bin_map) assert np.sum(np.isnan(new_uvw)) == 0 + assert np.sum(np.isnan(new_times)) == 0 - new_times = np.array([dct_inv[ii][2] for ii in range(len(atset))]) - new_times = np.mean(np.array(list_of_timebins), axis=1)[new_times] - + new_times = new_times[np.array([dct_inv[ii][2] for ii in range(len(atset))])] assert np.all(np.diff(new_times) >= 0) new_ant1 = np.array([dct_inv[ii][0] for ii in range(len(atset))]) -- GitLab From 303fb5fa38b629bef9d1f92349ec8fa5f1444f29 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 29 Oct 2021 12:23:35 +0200 Subject: [PATCH 50/89] Observation: flag baselines more verbose --- resolve/data/observation.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/resolve/data/observation.py b/resolve/data/observation.py index 7c41e868..cc4e25d4 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -589,7 +589,11 @@ class Observation(BaseObservation): ind = np.logical_and(ant1 == ant1_index, ant2 == ant2_index) wgt = self._weight.copy() wgt[:, ind] = 0. - print(f"INFO: Flag baseline {ant1_index}-{ant2_index}, {np.sum(ind)}/{self.nrow} rows flagged.") + antenna_names = self.auxiliary_table("ANTENNA")["STATION"] + ant1_name = antenna_names[ant1_index] + ant2_name = antenna_names[ant2_index] + if np.sum(ind) > 0: + print(f"INFO: Flag baseline {ant1_name}-{ant2_name}, {np.sum(ind)}/{self.nrow} rows flagged.") return Observation(self._antpos, self._vis, wgt, self._polarization, self._freq, self._auxiliary_tables) @property -- GitLab From 54d7fadc088f0f1168376f5ef4de204ddbfe68da Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 29 Oct 2021 12:34:47 +0200 Subject: [PATCH 51/89] Observation: Fixup --- resolve/data/observation.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/resolve/data/observation.py b/resolve/data/observation.py index cc4e25d4..dd9299c3 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -433,7 +433,7 @@ class Observation(BaseObservation): def __getitem__(self, slc, copy=False): # FIXME Do I need to change something in self._auxiliary_tables? ap = self._antpos[slc] - vis = self._vis[slc] + vis = self._vis[:, slc] wgt = self._weight[:, slc] if copy: ap = ap.copy() @@ -556,6 +556,8 @@ class Observation(BaseObservation): enum = np.bincount(masterindex, weights=vis[pol, :, freq].real*wgt[pol, :, freq]) enum = enum + 1j*np.bincount(masterindex, weights=vis[pol, :, freq].imag*wgt[pol, :, freq]) denom = np.bincount(masterindex, weights=wgt[pol, :, freq]) + if np.min(denom) == 0.: + raise ValueError("Time bin with total weight 0. detected.") new_vis[pol, :, freq] = enum/denom new_wgt[pol, :, freq] = denom -- GitLab From 9885801cba37b5f6239b0499c71b12e8c22dec8f Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 29 Oct 2021 16:23:11 +0200 Subject: [PATCH 52/89] Observation: restrict_to_stokes_i does nothing if no polarization is there --- resolve/data/observation.py | 1 - 1 file changed, 1 deletion(-) diff --git a/resolve/data/observation.py b/resolve/data/observation.py index dd9299c3..d3f13ea6 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -497,7 +497,6 @@ class Observation(BaseObservation): def restrict_to_stokesi(self): # FIXME Do I need to change something in self._auxiliary_tables? - my_asserteq(self._vis.shape[0], 4) ind = self._polarization.stokes_i_indices() vis = self._vis[ind] wgt = self._weight[ind] -- GitLab From a2b2d77e9ae6121202083c333319a191a41ebeed Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Mon, 8 Nov 2021 10:27:14 +0100 Subject: [PATCH 53/89] Upgrade demo to new nifty possibilities --- demo/cygnusa.cfg | 31 ++ demo/imaging_with_automatic_weighting.py | 358 ++++++++++++----------- resolve/simple_operators.py | 13 + test/test_operators.py | 21 ++ 4 files changed, 253 insertions(+), 170 deletions(-) create mode 100644 demo/cygnusa.cfg create mode 100644 test/test_operators.py diff --git a/demo/cygnusa.cfg b/demo/cygnusa.cfg new file mode 100644 index 00000000..d70753c4 --- /dev/null +++ b/demo/cygnusa.cfg @@ -0,0 +1,31 @@ +[data] +path = CYG-ALL-2052-2MHZ.ms + +[output] +directory = Cygnus_A_2052MHz + +[response] +wgridding = False +epsilon = 1e-6 + +[sky] +space npix x = 1000 +space npix y = 1000 +space fov x = 0.05deg +space fov y = 0.05deg + +zero mode offset = 21 +zero mode mean = 1 +zero mode stddev = 0.1 + +space fluctuations mean = 5 +space fluctuations stddev = 1 +space loglogavgslope mean = -2 +space loglogavgslope stddev = 0.5 +space flexibility mean = 1.2 +space flexibility stddev = 0.4 +space asperity mean = 0.2 +space asperity stddev = 0.2 + +[technical] +nthreads = 8 diff --git a/demo/imaging_with_automatic_weighting.py b/demo/imaging_with_automatic_weighting.py index e62a5c66..d4402e37 100644 --- a/demo/imaging_with_automatic_weighting.py +++ b/demo/imaging_with_automatic_weighting.py @@ -1,193 +1,211 @@ # SPDX-License-Identifier: GPL-3.0-or-later -# Copyright(C) 2019-2020 Max-Planck-Society -# Author: Philipp Arras - -import argparse +# Copyright(C) 2021 Max-Planck-Society +# Author: Philipp Arras, Jakob Knollmüller import numpy as np +import configparser +import os +import sys +from distutils.util import strtobool import nifty8 as ift import resolve as rve -from os.path import isfile, splitext - - -def main(): - parser = argparse.ArgumentParser() - # TODO OU process for automatic weighting - parser.add_argument("-j", type=int, default=1) - parser.add_argument("--use-cached", action="store_true") - parser.add_argument("--use-wgridding", action="store_true") - parser.add_argument( - "--data-column", - default="DATA", - help="Only active if a measurement set is read.", - ) - parser.add_argument("--point", action="append", nargs=2) - parser.add_argument("ms", type=str) - parser.add_argument("xfov", type=str) - parser.add_argument("yfov", type=str) - parser.add_argument("xpix", type=int) - parser.add_argument("ypix", type=int) - parser.add_argument("diffusefluxlevel", type=float) - args = parser.parse_args() - - ############################################################################ - # Define likelihood(s) - ############################################################################ - if splitext(args.ms)[1] == ".npz": - obs = rve.Observation.load(args.ms) - else: - obs = rve.ms2observations(args.ms, args.data_column, False, 0, "stokesiavg")[0] - - rve.set_nthreads(args.j) - rve.set_wgridding(args.use_wgridding) - fov = np.array([rve.str2rad(args.xfov), rve.str2rad(args.yfov)]) - npix = np.array([args.xpix, args.ypix]) - rve.set_epsilon(1 / 10 / obs.max_snr()) - - dom = ift.RGSpace(npix, fov / npix) - logsky = ift.SimpleCorrelatedField( - dom, args.diffusefluxlevel, (1, 0.1), (5, 1), (1.2, 0.4), (0.2, 0.2), (-2, 0.5) - ) - diffuse = logsky.exp() - if args.point is not None: - ppos = [] - for point in args.point: - ppos.append([rve.str2rad(point[0]), rve.str2rad(point[1])]) - inserter = rve.PointInserter(dom, ppos) - points = ift.InverseGammaOperator( - inserter.domain, alpha=0.5, q=0.2 / dom.scalar_dvol - ).ducktape("points") - points = inserter @ points - sky = diffuse + points +from matplotlib.colors import LogNorm +from resolve.mpi import master + + +def main(cfg_file_name): + cfg = configparser.ConfigParser() + cfg.read(cfg_file_name) + if master: + print(f"Load {cfg_file_name}") + + point_sources = [("0deg", "0deg")] + + rve.set_epsilon(cfg["response"].getfloat("epsilon")) + rve.set_wgridding(strtobool(cfg["response"]["wgridding"])) + rve.set_double_precision(True) + rve.set_nthreads(cfg["technical"].getint("nthreads")) + + nx = cfg["sky"].getint("space npix x") + ny = cfg["sky"].getint("space npix y") + dx = rve.str2rad(cfg["sky"]["space fov x"]) / nx + dy = rve.str2rad(cfg["sky"]["space fov y"]) / ny + dom = ift.RGSpace([nx, ny], [dx, dy]) + + # Sky model + # Diffuse + cfm_zm_args = { + "offset_mean": + cfg["sky"].getfloat("zero mode offset"), + "offset_std": ( + cfg["sky"].getfloat("zero mode mean"), + cfg["sky"].getfloat("zero mode stddev"), + ), + } + cfm_spatial_args = { + "fluctuations": ( + cfg["sky"].getfloat("space fluctuations mean"), + cfg["sky"].getfloat("space fluctuations stddev"), + ), + "loglogavgslope": ( + cfg["sky"].getfloat("space loglogavgslope mean"), + cfg["sky"].getfloat("space loglogavgslope stddev"), + ), + "flexibility": ( + cfg["sky"].getfloat("space flexibility mean"), + cfg["sky"].getfloat("space flexibility stddev"), + ), + "asperity": ( + cfg["sky"].getfloat("space asperity mean"), + cfg["sky"].getfloat("space asperity stddev"), + ), + } + logdiffuse = ift.SimpleCorrelatedField(dom, **cfm_spatial_args, **cfm_zm_args) + + # Point sources + ppos = [] + for point in point_sources: + ppos.append([rve.str2rad(point[0]), rve.str2rad(point[1])]) + inserter = rve.PointInserter(dom, ppos) + points = ift.InverseGammaOperator( + inserter.domain, alpha=0.5, q=0.2 / dom.scalar_dvol).ducktape("points") + points = inserter @ points + + sky = logdiffuse.exp() + points + # /Sky model + + p = ift.Plot() + for _ in range(9): + p.add(sky(ift.from_random(sky.domain)), norm=LogNorm()) + if master: + p.output(name="sky_prior_samples.png") + + obs_file_name = cfg["data"]["path"] + if os.path.splitext(obs_file_name)[1] == ".npz": + obs = rve.Observation.load(obs_file_name).restrict_to_stokesi().average_stokesi() + elif os.path.splitext(obs_file_name)[1] == ".ms": + obs = rve.ms2observations(obs_file_name, "DATA", False, 0, "stokesiavg") + assert len(obs) == 1 + obs = obs[0] else: - sky = diffuse - # TODO Add mode with independent noise learning + raise RuntimeError( + f"Do not understand file name ending of {obs_file_name}. Supported: ms, npz.") + + # Weightop npix = 2500 effuv = obs.effective_uvwlen().val[0] assert obs.nfreq == obs.npol == 1 dom = ift.RGSpace(npix, 2 * np.max(effuv) / npix) - logwgt = ift.SimpleCorrelatedField( - dom, 0, (2, 2), (2, 2), (1.2, 0.4), (0.5, 0.2), (-2, 0.5), "invcov" - ) - li = ift.LinearInterpolator(dom, effuv) - weightop = ift.makeOp(obs.weight) @ ( - rve.AddEmptyDimension(li.target) @ li @ logwgt.exp() - ) ** (-2) - - plotter = rve.Plotter("png", "plots") - plotter.add("logsky", logsky) - plotter.add("power spectrum logsky", logsky.power_spectrum) - plotter.add("bayesian weighting", logwgt.exp()) - plotter.add("power spectrum bayesian weighting", logwgt.power_spectrum) - - ############################################################################ - # MINIMIZATION - ############################################################################ - if rve.mpi.master: - if args.point is not None: - # MAP points with original weights - lh = rve.ImagingLikelihood(obs, points) - ham = ift.StandardHamiltonian(lh) - state = rve.MinimizationState(0.1 * ift.from_random(ham.domain), []) - mini = ift.NewtonCG( - ift.GradientNormController(name="newton", iteration_limit=4) - ) - if args.use_cached and isfile("stage0"): - state = rve.MinimizationState.load("stage0") - else: - state = rve.simple_minimize(ham, state.mean, 0, mini) - plotter.plot("stage0", state) - state.save("stage0") - - # MAP diffuse with original weights + logwgt = ift.SimpleCorrelatedField(dom, 0, (2, 2), (2, 2), (1.2, 0.4), (0.5, 0.2), (-2, 0.5), + "invcov") + li = ift.LinearInterpolator(dom, effuv.T) + conv = rve.DomainChangerAndReshaper(li.target, obs.weight.domain) + weightop = ift.makeOp(obs.weight) @ (conv @ li @ logwgt.exp()) ** (-2) + # /Weightop + + full_lh = rve.ImagingLikelihood(obs, sky, inverse_covariance_operator=weightop) + position = 0.1 * ift.from_random(full_lh.domain) + + common = { + "output_directory": cfg["output"]["directory"], + "plottable_operators": { + "sky": sky, + "logsky": sky.log(), + "points": points, + "logdiffuse": logdiffuse, + "logdiffuse pspec": logdiffuse.power_spectrum + }, + "overwrite": True + } + + if master: + # Points only, MAP + lh = rve.ImagingLikelihood(obs, points) + sl = ift.optimize_kl( + lh, + 1, + 0, + ift.NewtonCG(ift.GradientNormController(name="hamiltonian", iteration_limit=4)), + None, + None, + initial_position=position, + **common) + position = ift.MultiField.union([position, sl.local_item(0)]) + + # Points constant, diffuse only, MAP lh = rve.ImagingLikelihood(obs, sky) - plotter.add_histogram( - "normalized residuals (original weights)", lh.normalized_residual - ) - ham = ift.StandardHamiltonian(lh) - if args.point is None: - fld = 0.1 * ift.from_random(diffuse.domain) - else: - fld = ift.MultiField.union( - [0.1 * ift.from_random(diffuse.domain), state.mean] - ) - state = rve.MinimizationState(fld, []) - mini = ift.NewtonCG( - ift.GradientNormController(name="newton", iteration_limit=20) - ) - if args.use_cached and isfile("stage1"): - state = rve.MinimizationState.load("stage1") - else: - state = rve.simple_minimize(ham, state.mean, 0, mini) - plotter.plot("stage1", state) - state.save("stage1") - - # Only weights - lh = rve.ImagingLikelihood(obs, sky, inverse_covariance_operator=weightop) - plotter.add_histogram( - "normalized residuals (learned weights)", lh.normalized_residual - ) - ic = ift.AbsDeltaEnergyController(0.1, 3, 100, name="Sampling") - ham = ift.StandardHamiltonian(lh, ic) + sl = ift.optimize_kl( + lh, + 1, + 0, + ift.NewtonCG(ift.GradientNormController(name="hamiltonian", iteration_limit=20)), + None, + None, + initial_position=position, + initial_index=1, + **common) + position = ift.MultiField.union([position, sl.local_item(0)]) + cst = sky.domain.keys() - state = rve.MinimizationState( - ift.MultiField.union([0.1 * ift.from_random(weightop.domain), state.mean]), - [], - ) - mini = ift.VL_BFGS(ift.GradientNormController(name="bfgs", iteration_limit=20)) - if args.use_cached and isfile("stage2"): - state = rve.MinimizationState.load("stage2") - else: - for ii in range(10): - state = rve.simple_minimize(ham, state.mean, 0, mini, cst, cst) - plotter.plot(f"stage2_{ii}", state) - state.save("stage2") + sl = ift.optimize_kl( + full_lh, + 1, + 0, + ift.VL_BFGS(ift.GradientNormController(name="bfgs", iteration_limit=20)), + ift.AbsDeltaEnergyController(0.1, 3, 100, name="Sampling"), + None, + constants=cst, + point_estimates=cst, + initial_position=position, + initial_index=2, + **common) + position = ift.MultiField.union([position, sl.local_item(0)]) if rve.mpi.mpi: - if not rve.mpi.master: - state = None - state = rve.mpi.comm.bcast(state, root=0) + if not master: + position = None + position = rve.mpi.comm.bcast(position, root=0) ift.random.push_sseq_from_seed(42) - # MGVI sky - ic = ift.AbsDeltaEnergyController(0.1, 3, 200, name="Sampling") - lh = rve.ImagingLikelihoodVariableCovariance(obs, sky, weightop) - ham = ift.StandardHamiltonian(lh, ic) - if args.point is not None: - cst = list(points.domain.keys()) + list(weightop.domain.keys()) - else: - cst = list(weightop.domain.keys()) - mini = ift.NewtonCG(ift.GradientNormController(name="newton", iteration_limit=15)) - for ii in range(4): - fname = f"stage3_{ii}" - if args.use_cached and isfile(fname): - state = rve.MinimizationState.load(fname) + # First MGVI diffuse sky, then sky + weighting simultaneously + + def get_mini(iglobal): + if iglobal < 7: + return ift.NewtonCG(ift.GradientNormController(name="kl", iteration_limit=15)) + elif iglobal < 7 + 5: + return ift.VL_BFGS(ift.GradientNormController(name="kl", iteration_limit=15)) else: - state = rve.simple_minimize(ham, state.mean, 5, mini, cst, cst) - plotter.plot(f"stage3_{ii}", state) - state.save(fname) - - # Sky + weighting simultaneously - ic = ift.AbsDeltaEnergyController(0.1, 3, 700, name="Sampling") - ham = ift.StandardHamiltonian(lh, ic) - for ii in range(30): - if ii < 5: - mini = ift.VL_BFGS( - ift.GradientNormController(name="newton", iteration_limit=15) - ) + return ift.NewtonCG(ift.GradientNormController(name="kl", iteration_limit=15)) + + def get_sampling(iglobal): + if iglobal < 7: + lim = 200 + elif iglobal < 7 + 5: + lim = 700 else: - mini = ift.NewtonCG( - ift.GradientNormController(name="newton", iteration_limit=15) - ) - fname = f"stage4_{ii}" - if args.use_cached and isfile(fname): - state = rve.MinimizationState.load(fname) + lim = 700 + return ift.AbsDeltaEnergyController(deltaE=0.5, iteration_limit=lim) + + def get_cst(iglobal): + if iglobal < 7: + return list(points.domain.keys()) + list(weightop.domain.keys()) else: - state = rve.simple_minimize(ham, state.mean, 5, mini) - plotter.plot(f"stage4_{ii}", state) - state.save(fname) + return [] + + sl = ift.optimize_kl( + full_lh, + 35, + 5, + get_mini, + get_sampling, + None, + constants=get_cst, + point_estimates=get_cst, + initial_position=position, + initial_index=7, + **common) if __name__ == "__main__": - main() + main(sys.argv[1]) diff --git a/resolve/simple_operators.py b/resolve/simple_operators.py index 44afc895..8be64027 100644 --- a/resolve/simple_operators.py +++ b/resolve/simple_operators.py @@ -113,3 +113,16 @@ def MultiDomainVariableCovarianceGaussianEnergy(data, signal_response, invcov): resi = KeyPrefixer(data.domain, "resi") @ ift.Adder(data, True) @ signal_response invcov = KeyPrefixer(data.domain, "icov") @ invcov return reduce(add, res) @ (resi + invcov) + + +class DomainChangerAndReshaper(ift.LinearOperator): + def __init__(self, domain, target): + self._domain = ift.DomainTuple.make(domain) + self._target = ift.DomainTuple.make(target) + self._capability = self.TIMES | self.ADJOINT_TIMES + + def apply(self, x, mode): + self._check_input(x, mode) + x = x.val + tgt = self._tgt(mode) + return ift.makeField(tgt, x.reshape(tgt.shape)) diff --git a/test/test_operators.py b/test/test_operators.py new file mode 100644 index 00000000..b4c9daa9 --- /dev/null +++ b/test/test_operators.py @@ -0,0 +1,21 @@ +# SPDX-License-Identifier: GPL-3.0-or-later +# Copyright(C) 2020-2021 Max-Planck-Society +# Author: Philipp Arras + +from os.path import join + +import numpy as np +import pytest + +import nifty8 as ift +import resolve as rve + +pmp = pytest.mark.parametrize +np.seterr(all="raise") + + +def test_reshaper(): + dom = ift.UnstructuredDomain([3]), ift.RGSpace([5]) + tgt = ift.RGSpace([15]) + op = rve.DomainChangerAndReshaper(dom, tgt) + ift.extra.check_linear_operator(op) -- GitLab From 8cccee0354d22041d219dce1f9cb4672bc4684ef Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Tue, 9 Nov 2021 13:35:13 +0100 Subject: [PATCH 54/89] Add weighting and point sources to config file --- demo/cygnusa.cfg | 24 +- demo/imaging_with_automatic_weighting.py | 212 +++++-------- resolve/__init__.py | 3 +- resolve/minimization.py | 137 --------- resolve/mpi.py | 2 + resolve/plotter.py | 375 ----------------------- resolve/sky_model.py | 67 ++++ 7 files changed, 170 insertions(+), 650 deletions(-) delete mode 100644 resolve/minimization.py delete mode 100644 resolve/plotter.py create mode 100644 resolve/sky_model.py diff --git a/demo/cygnusa.cfg b/demo/cygnusa.cfg index d70753c4..c9f8bb57 100644 --- a/demo/cygnusa.cfg +++ b/demo/cygnusa.cfg @@ -27,5 +27,27 @@ space flexibility stddev = 0.4 space asperity mean = 0.2 space asperity stddev = 0.2 +point sources mode = single +point sources locations = 0deg,0deg;0.35as,-0.22as +point sources alpha = 0.5 +point sources q = 0.2 + +[weighting] +enable = True +model = cfm + +npix = 2500 +zero mode offset = 0 +zero mode mean = 2 +zero mode stddev = 2 +fluctuations mean = 2 +fluctuations stddev = 2 +loglogavgslope mean = -2 +loglogavgslope stddev = 0.5 +flexibility mean = 1.2 +flexibility stddev = 0.4 +asperity mean = 0.5 +asperity stddev = 0.2 + [technical] -nthreads = 8 +nthreads = 4 diff --git a/demo/imaging_with_automatic_weighting.py b/demo/imaging_with_automatic_weighting.py index d4402e37..cf0f3038 100644 --- a/demo/imaging_with_automatic_weighting.py +++ b/demo/imaging_with_automatic_weighting.py @@ -1,14 +1,14 @@ # SPDX-License-Identifier: GPL-3.0-or-later # Copyright(C) 2021 Max-Planck-Society -# Author: Philipp Arras, Jakob Knollmüller +# Author: Philipp Arras -import numpy as np import configparser import os import sys from distutils.util import strtobool import nifty8 as ift +import numpy as np import resolve as rve from matplotlib.colors import LogNorm from resolve.mpi import master @@ -20,67 +20,14 @@ def main(cfg_file_name): if master: print(f"Load {cfg_file_name}") - point_sources = [("0deg", "0deg")] - rve.set_epsilon(cfg["response"].getfloat("epsilon")) rve.set_wgridding(strtobool(cfg["response"]["wgridding"])) rve.set_double_precision(True) rve.set_nthreads(cfg["technical"].getint("nthreads")) - nx = cfg["sky"].getint("space npix x") - ny = cfg["sky"].getint("space npix y") - dx = rve.str2rad(cfg["sky"]["space fov x"]) / nx - dy = rve.str2rad(cfg["sky"]["space fov y"]) / ny - dom = ift.RGSpace([nx, ny], [dx, dy]) - - # Sky model - # Diffuse - cfm_zm_args = { - "offset_mean": - cfg["sky"].getfloat("zero mode offset"), - "offset_std": ( - cfg["sky"].getfloat("zero mode mean"), - cfg["sky"].getfloat("zero mode stddev"), - ), - } - cfm_spatial_args = { - "fluctuations": ( - cfg["sky"].getfloat("space fluctuations mean"), - cfg["sky"].getfloat("space fluctuations stddev"), - ), - "loglogavgslope": ( - cfg["sky"].getfloat("space loglogavgslope mean"), - cfg["sky"].getfloat("space loglogavgslope stddev"), - ), - "flexibility": ( - cfg["sky"].getfloat("space flexibility mean"), - cfg["sky"].getfloat("space flexibility stddev"), - ), - "asperity": ( - cfg["sky"].getfloat("space asperity mean"), - cfg["sky"].getfloat("space asperity stddev"), - ), - } - logdiffuse = ift.SimpleCorrelatedField(dom, **cfm_spatial_args, **cfm_zm_args) - - # Point sources - ppos = [] - for point in point_sources: - ppos.append([rve.str2rad(point[0]), rve.str2rad(point[1])]) - inserter = rve.PointInserter(dom, ppos) - points = ift.InverseGammaOperator( - inserter.domain, alpha=0.5, q=0.2 / dom.scalar_dvol).ducktape("points") - points = inserter @ points - - sky = logdiffuse.exp() + points - # /Sky model - - p = ift.Plot() - for _ in range(9): - p.add(sky(ift.from_random(sky.domain)), norm=LogNorm()) - if master: - p.output(name="sky_prior_samples.png") + enable_weighting = cfg["weighting"].getboolean("enable") + # Data obs_file_name = cfg["data"]["path"] if os.path.splitext(obs_file_name)[1] == ".npz": obs = rve.Observation.load(obs_file_name).restrict_to_stokesi().average_stokesi() @@ -91,84 +38,85 @@ def main(cfg_file_name): else: raise RuntimeError( f"Do not understand file name ending of {obs_file_name}. Supported: ms, npz.") + # /Data + + # Sky model + sky, operators = rve.single_frequency_sky(cfg["sky"]) + enable_points = operators["points"] is not None + operators["sky"] = sky + operators["logsky"] = sky.log() + + p = ift.Plot() + for _ in range(9): + p.add(sky(ift.from_random(sky.domain)), norm=LogNorm()) + if master: + p.output(name="sky_prior_samples.png") + # /Sky model - # Weightop - npix = 2500 - effuv = obs.effective_uvwlen().val[0] - assert obs.nfreq == obs.npol == 1 - dom = ift.RGSpace(npix, 2 * np.max(effuv) / npix) - logwgt = ift.SimpleCorrelatedField(dom, 0, (2, 2), (2, 2), (1.2, 0.4), (0.5, 0.2), (-2, 0.5), - "invcov") - li = ift.LinearInterpolator(dom, effuv.T) - conv = rve.DomainChangerAndReshaper(li.target, obs.weight.domain) - weightop = ift.makeOp(obs.weight) @ (conv @ li @ logwgt.exp()) ** (-2) - # /Weightop + # Bayesian weighting + if enable_weighting: + assert obs.nfreq == obs.npol == 1 + subcfg = cfg["weighting"] + if subcfg["model"] == "cfm": + npix = 2500 + xs = np.log(obs.effective_uvwlen().val) + xs -= np.min(xs) + if not xs.shape[0] == xs.shape[2] == 1: + raise RuntimeError + # FIXME Use Integrated Wiener process + dom = ift.RGSpace(npix, 2 * np.max(xs) / npix) + logwgt, cfm = rve.cfm_from_cfg(subcfg, dom, "invcov") + li = ift.LinearInterpolator(dom, xs[0].T) + conv = rve.DomainChangerAndReshaper(li.target, obs.weight.domain) + weightop = ift.makeOp(obs.weight) @ (conv @ li @ logwgt.exp()) ** (-2) + operators["logweights"] = logwgt + operators["weights"] = logwgt.exp() + operators["logweights power spectrum"] = cfm.power_spectrum + else: + raise NotImplementedError + else: + weightop = None + # /Bayesian weighting full_lh = rve.ImagingLikelihood(obs, sky, inverse_covariance_operator=weightop) position = 0.1 * ift.from_random(full_lh.domain) - - common = { - "output_directory": cfg["output"]["directory"], - "plottable_operators": { - "sky": sky, - "logsky": sky.log(), - "points": points, - "logdiffuse": logdiffuse, - "logdiffuse pspec": logdiffuse.power_spectrum - }, - "overwrite": True - } + common = {"output_directory": cfg["output"]["directory"], + "plottable_operators": operators, + "overwrite": True} if master: - # Points only, MAP - lh = rve.ImagingLikelihood(obs, points) - sl = ift.optimize_kl( - lh, - 1, - 0, - ift.NewtonCG(ift.GradientNormController(name="hamiltonian", iteration_limit=4)), - None, - None, - initial_position=position, - **common) - position = ift.MultiField.union([position, sl.local_item(0)]) + common["comm"] = rve.mpi.comm_self + if enable_points: + # Points only, MAP + lh = rve.ImagingLikelihood(obs, operators["points"]) + mini = ift.NewtonCG(ift.GradientNormController(name="hamiltonian", iteration_limit=4)) + sl = ift.optimize_kl(lh, 1, 0, mini, None, None, initial_position=position, **common) + position = sl.local_item(0) # Points constant, diffuse only, MAP lh = rve.ImagingLikelihood(obs, sky) - sl = ift.optimize_kl( - lh, - 1, - 0, - ift.NewtonCG(ift.GradientNormController(name="hamiltonian", iteration_limit=20)), - None, - None, - initial_position=position, - initial_index=1, - **common) - position = ift.MultiField.union([position, sl.local_item(0)]) + mini = ift.NewtonCG(ift.GradientNormController(name="hamiltonian", iteration_limit=20)) + sl = ift.optimize_kl(lh, 1, 0, mini, None, None, initial_position=position, initial_index=1, **common) + position = sl.local_item(0) + # Only weighting, MGVI cst = sky.domain.keys() - sl = ift.optimize_kl( - full_lh, - 1, - 0, - ift.VL_BFGS(ift.GradientNormController(name="bfgs", iteration_limit=20)), - ift.AbsDeltaEnergyController(0.1, 3, 100, name="Sampling"), - None, - constants=cst, - point_estimates=cst, - initial_position=position, - initial_index=2, - **common) - position = ift.MultiField.union([position, sl.local_item(0)]) + mini = ift.VL_BFGS(ift.GradientNormController(name="bfgs", iteration_limit=20)) + ic_sampling = ift.AbsDeltaEnergyController(0.1, 3, 100, name="Sampling") + sl = ift.optimize_kl(full_lh, 5, 3, mini, ic_sampling, None, + constants=cst, point_estimates=cst, + initial_position=position, initial_index=2, **common) + position = sl.local_item(0) + + # Reset sky + position = ift.MultiField.union([0.1*ift.from_random(sky.domain), position.extract(weightop.domain)]) if rve.mpi.mpi: if not master: position = None position = rve.mpi.comm.bcast(position, root=0) ift.random.push_sseq_from_seed(42) - - # First MGVI diffuse sky, then sky + weighting simultaneously + common["comm"] = rve.mpi.comm def get_mini(iglobal): if iglobal < 7: @@ -185,26 +133,20 @@ def main(cfg_file_name): lim = 700 else: lim = 700 - return ift.AbsDeltaEnergyController(deltaE=0.5, iteration_limit=lim) + return ift.AbsDeltaEnergyController(deltaE=0.5, iteration_limit=lim, name="Sampling") def get_cst(iglobal): + res = [] if iglobal < 7: - return list(points.domain.keys()) + list(weightop.domain.keys()) - else: - return [] - - sl = ift.optimize_kl( - full_lh, - 35, - 5, - get_mini, - get_sampling, - None, - constants=get_cst, - point_estimates=get_cst, - initial_position=position, - initial_index=7, - **common) + res += list(weightop.domain.keys()) + if enable_points: + res += list(operators["points"].domain.keys()) + return res + + # First MGVI diffuse sky, then sky + weighting simultaneously + sl = ift.optimize_kl(full_lh, 35, 6, get_mini, get_sampling, None, + constants=get_cst, point_estimates=get_cst, + initial_position=position, initial_index=7, **common) if __name__ == "__main__": diff --git a/resolve/__init__.py b/resolve/__init__.py index 58bb5754..8321f28e 100644 --- a/resolve/__init__.py +++ b/resolve/__init__.py @@ -8,7 +8,6 @@ from .data.averaging import * from .fits import field2fits, fits2field from .global_config import * from .likelihood import * -from .minimization import Minimization, MinimizationState, simple_minimize from .mosaicing import * from .mpi import onlymaster from .mpi_operators import * @@ -20,7 +19,6 @@ from .multi_frequency.operators import ( WienerIntegrations, ) from .data.observation import * -from .plotter import MfPlotter, Plotter from .points import PointInserter from .data.polarization import Polarization from .polarization_matrix_exponential import * @@ -29,3 +27,4 @@ from .response import MfResponse, ResponseDistributor, StokesIResponse, SingleRe from .simple_operators import * from .util import * from .extra import mpi_load +from .sky_model import * diff --git a/resolve/minimization.py b/resolve/minimization.py deleted file mode 100644 index dbf004c5..00000000 --- a/resolve/minimization.py +++ /dev/null @@ -1,137 +0,0 @@ -# SPDX-License-Identifier: GPL-3.0-or-later -# Copyright(C) 2020-2021 Max-Planck-Society -# Author: Philipp Arras - -import pickle - -import nifty8 as ift - -from .mpi import comm, onlymaster -from .util import compare_attributes, my_assert, my_assert_isinstance, my_asserteq - - -def simple_minimize( - operator, position, n_samples, minimizer, constants=[], point_estimates=[] -): - mini = Minimization(operator, position, n_samples, constants, point_estimates, comm) - return mini.minimize(minimizer) - - -class Minimization: - def __init__(self, operator, position, n_samples, constants=[], point_estimates=[], comm=None): - n_samples = int(n_samples) - self._position = position - position = position.extract(operator.domain) - if n_samples == 0: - self._e = ift.EnergyAdapter(position, operator, constants, True, True) - self._n = n_samples - else: - my_assert(n_samples > 0) - dct = { - "position": position, - "hamiltonian": operator, - "n_samples": n_samples, - "minimizer_sampling": None, - "constants": constants, - "point_estimates": point_estimates, - "mirror_samples": True, - "comm": comm, - "nanisinf": True, - } - self._e = ift.SampledKLEnergy(**dct) - self._n, self._m = dct["n_samples"], dct["mirror_samples"] - - def minimize(self, minimizer): - self._e, _ = minimizer(self._e) - position = ift.MultiField.union([self._position, self._e.position]) - if self._n == 0: - return MinimizationState(position, []) - zeros = ift.full(position.domain, 0.) - samples = [ss.unite(zeros) for ss in self._e.samples] - my_asserteq(len(samples), 2 * self._n if self._m else self._n) - return MinimizationState(position, samples, self._m) - - -class MinimizationState: - def __init__(self, position, samples=[], mirror_samples=False): - self._samples = list(samples) - self._position = position - if len(samples) > 0: - my_asserteq(samples[0].domain, *[ss.domain for ss in samples]) - my_assert(set(samples[0].domain.keys()) <= set(position.domain.keys())) - self._mirror = bool(mirror_samples) - - @property - def has_samples(self): - return len(self._samples) > 0 - - def __getitem__(self, key): - # FIXME Add MPI support - if not isinstance(key, int): - raise TypeError - if key >= len(self) or key < 0: - raise IndexError - if key == 0 and not self.has_samples: - return self._position - if self._mirror and key >= len(self) // 2: - return self._position.unite(-self._samples[key]) - return self._position.unite(self._samples[key]) - - def all_samples(self): - if not self.has_samples: - return None - lst = [] - if self._mirror: - lst = [-ss for ss in self._samples] - return lst + self._samples - - def __len__(self): - if not self.has_samples: - return 1 - return (2 if self._mirror else 1) * len(self._samples) - - def __eq__(self, other): - if not isinstance(other, MinimizationState): - return False - return compare_attributes(self, other, ("_samples", "_position", "_mirror")) - - @property - def mean(self): - return self._position - - @property - def domain(self): - return self._position.domain - - @onlymaster - def save(self, file_name): - with open(file_name, "wb") as f: - pickle.dump( - [self._position, self._samples, self._mirror], - f, - pickle.HIGHEST_PROTOCOL, - ) - - @staticmethod - def load(file_name): - with open(file_name, "rb") as f: - position, samples, mirror = pickle.load(f) - my_assert_isinstance(position, (ift.MultiField, ift.Field)) - my_assert_isinstance(mirror, bool) - return MinimizationState(position, samples, mirror) - - - def operator_stats(self, operator): - """Return ift.StatCalculator for operator applied to all samples""" - sc = ift.StatCalculator() - for ss in self: - sc.add(operator.force(ss)) - return sc - - - def operator_samples(self, operator): - """Return posterior samples for operator""" - res = [] - for ss in self: - res.append(operator.force(ss)) - return res diff --git a/resolve/mpi.py b/resolve/mpi.py index 73507781..0328336f 100644 --- a/resolve/mpi.py +++ b/resolve/mpi.py @@ -11,6 +11,7 @@ try: master = MPI.COMM_WORLD.Get_rank() == 0 comm = MPI.COMM_WORLD + comm_self = MPI.COMM_SELF ntask = comm.Get_size() rank = comm.Get_rank() master = rank == 0 @@ -19,6 +20,7 @@ except ImportError: master = True mpi = False comm = None + comm_self = None rank = 0 diff --git a/resolve/plotter.py b/resolve/plotter.py deleted file mode 100644 index 90c09dcc..00000000 --- a/resolve/plotter.py +++ /dev/null @@ -1,375 +0,0 @@ -# SPDX-License-Identifier: GPL-3.0-or-later -# Copyright(C) 2020 Max-Planck-Society -# Author: Philipp Arras - -import os -from os.path import join, splitext - -import matplotlib.pyplot as plt -import numpy as np - -import nifty8 as ift - -from .minimization import MinimizationState -from .mpi import onlymaster -from .util import my_assert_isinstance, my_asserteq - -UNIT = 6 - - -@onlymaster -def makedirs(*args, **kwargs): - from os import makedirs - - makedirs(*args, **kwargs) - - -class Plotter: - def __init__(self, fileformat, directory): - self._nifty, self._calib, self._hist = [], [], [] - self._fits = [] - self._f = fileformat - self._dir = directory - makedirs(self._dir, exist_ok=True) - - @onlymaster - def add_fits(self, name, operator, direction=None): - self._fits.append({"operator": operator, "title": str(name), - "direction": direction}) - - @onlymaster - def add_npy(self, name, operator, nsamples=0): - """Save output of operator as npy files. - - This function may be useful for saving results of a resolve run and - plotting them later with dedicated plotting scripts. - - If the plotter gets samples during plotting, the posterior mean and the - posterior standard deviation are saved. Additionally, a number of - posterior samples are saved. If no samples are available, only one npy - array will be saved. - - Parameters - ---------- - name : str - Naming of the output directory and files. - operator : Operator - The output of this operator is saved - nsamples : int - Determines how many posterior samples shall be saved additionally to - the posterior mean. - - """ - raise NotImplementedError - - @onlymaster - def add(self, name, operator, **kwargs): - my_assert_isinstance(operator, ift.Operator) - self._nifty.append({"operator": operator, "title": str(name), "kwargs": kwargs}) - - @onlymaster - def add_calibration_solution(self, name, operator): - my_assert_isinstance(operator, ift.Operator) - self._calib.append({"title": str(name), "operator": operator}) - - @onlymaster - def add_histogram(self, name, operator): - my_assert_isinstance(operator, ift.Operator) - self._hist.append({"title": str(name), "operator": operator}) - - @onlymaster - def plot(self, identifier, state): - my_assert_isinstance(state, (ift.MultiField, MinimizationState)) - for ii, obj in enumerate(self._nifty): - op, fname = self._plot_init(obj, state, identifier) - if fname is None: - continue - try: - _plot_nifty(state, op, obj["kwargs"], fname) - except: - pass - for ii, obj in enumerate(self._calib): - op, fname = self._plot_init(obj, state, identifier) - if fname is None: - continue - _plot_calibration(state, op, fname) - for ii, obj in enumerate(self._hist): - op, fname = self._plot_init(obj, state, identifier) - if fname is None: - continue - _plot_histograms(state, fname, 20, 0.5, postop=op) - for ii, obj in enumerate(self._fits): - from .fits import field2fits - op, fname = self._plot_init(obj, state, identifier) - sc = ift.StatCalculator() - for ss in state: - sc.add(op.force(ss)) - field2fits(sc.mean, f"{fname}.fits", True, obj["direction"]) - mydir = join(self._dir, "zzz_latent") - makedirs(mydir, exist_ok=True) - _plot_histograms(state, join(mydir, f"{identifier}.{self._f}"), 5, 0.5) - - @onlymaster - def _plot_init(self, obj, state, identifier): - op = obj["operator"] - if not set(op.domain.keys()) <= set(state.domain.keys()): - return None, None - direc = join(self._dir, obj["title"]) - makedirs(direc, exist_ok=True) - fname = join(direc, f"{identifier}.{self._f}") - return op, fname - - @property - def directory(self): - return self._dir - - -class MfPlotter(Plotter): - def __init__(self, fileformat, directory): - super(MfPlotter, self).__init__(fileformat, directory) - self._spectra = [] - self._mf = [] - self._multi1d = [] - - @onlymaster - def add_spectra(self, name, operator, directions): - """ - Parameters - ---------- - directions: array, shape (n, 2) - In units of the space. - """ - my_assert_isinstance(operator, ift.Operator) - # FIXME Write interface checks - sdom = operator.target[1] - directions = np.round(np.array(directions) / np.array(sdom.distances)).astype( - int - ) - self._spectra.append( - {"title": str(name), "operator": operator, "directions": directions} - ) - - @onlymaster - def add_multiple2d(self, name, operator, directions=None, movie_length=2): - """ - Parameters - ---------- - directions: array, shape (n, 2) - In units of the space. - """ - my_assert_isinstance(operator, ift.Operator) - # FIXME Write interface checks - self._mf.append( - { - "title": str(name), - "operator": operator, - "directions": directions, - "movie_length": movie_length, - } - ) - - @onlymaster - def add_multiple1d(self, name, operator, **kwargs): - my_assert_isinstance(operator, ift.Operator) - my_asserteq(len(operator.target.shape), 2) - self._multi1d.append( - {"operator": operator, "title": str(name), "kwargs": kwargs} - ) - - @onlymaster - def plot(self, identifier, state): - super(MfPlotter, self).plot(identifier, state) - - for ii, obj in enumerate(self._spectra): - op, fname = self._plot_init(obj, state, identifier) - if fname is None: - continue - _plot_spectra(state, op, fname, obj["directions"]) - - for ii, obj in enumerate(self._mf): - op, fname = self._plot_init(obj, state, identifier) - if fname is None: - continue - _plot_mf(state, op, fname, obj["directions"], obj["movie_length"]) - - for ii, obj in enumerate(self._multi1d): - op, fname = self._plot_init(obj, state, identifier) - if fname is None: - continue - _plot_multi_oned(state, op, obj["kwargs"], fname) - - -def _plot_nifty(state, op, kwargs, fname): - if isinstance(state, MinimizationState) and len(state) > 0: - tgt = op.target - if ( - len(tgt) == 1 - and isinstance(tgt[0], (ift.RGSpace, ift.PowerSpace)) - and len(tgt.shape) == 1 - ): - p = ift.Plot() - p.add([op.force(ss) for ss in state], **kwargs) - p.output(xsize=UNIT, ysize=UNIT, name=fname) - else: - sc = ift.StatCalculator() - for ss in state: - sc.add(op.force(ss)) - p = ift.Plot() - p.add(sc.mean, **kwargs) - if len(state) > 1: - p.add(sc.var.sqrt() / sc.mean) - p.output(nx=2, ny=1, xsize=2 * UNIT, ysize=UNIT, name=fname) - else: - pos = state if isinstance(state, ift.MultiField) else state.mean - ift.single_plot(op.force(pos), **kwargs, name=fname) - - -def _plot_multi_oned(state, op, kwargs, fname): - if isinstance(state, MinimizationState) and len(state) > 0: - raise NotImplementedError - p = ift.Plot() - p.add([op.force(ss) for ss in state], **kwargs) - p.output(xsize=UNIT, ysize=UNIT, name=fname) - else: - pos = state if isinstance(state, ift.MultiField) else state.mean - ift.single_plot(op.force(pos), **kwargs, name=fname) - - -def _plot_calibration(state, op, fname): - npol, nants, _, nfreq = op.target.shape - if nfreq != 1: - raise NotImplementedError - lst = [op.force(ll) for ll in _generalstate2list(state)] - xs = np.arange(op.target[2].shape[0]) - if isinstance(op.target[2], ift.RGSpace): - xs = xs * op.target[2].distances[0] / 3600 - fig, axs = plt.subplots(2, 1, sharex=True, figsize=(4 * UNIT, npol * UNIT)) - axs = list(axs.ravel()) - for ii in range(npol): - axx = axs[ii] - axx.set_title(f"Polarization {ii}") - colors = plt.cm.viridis(np.linspace(0, 1, nants)) - for jj in range(nants): - for ll in lst: - axx.plot(xs, ll.val[ii, jj], alpha=0.3, color=colors[jj]) - axs[1].set_xlabel("Time [h]") - fig.savefig(fname) - plt.close(fig) - - -def _generalstate2list(state): - if isinstance(state, (ift.Field, ift.MultiField)): - return [state] - elif isinstance(state, MinimizationState): - if len(state) == 0: - return [state.mean] - else: - return state - raise TypeError - - -def _plot_histograms(state, fname, lim, width, postop=None): - lst = _generalstate2list(state) - if postop is not None: - lst0 = postop.force(lst[0]) - lst = (postop.force(ff) for ff in lst) - dom = postop.target - else: - lst0 = lst[0] - dom = state.domain - - if isinstance(dom, ift.DomainTuple): - iscplx = {"": np.iscomplexobj(lst0.val)} - else: - iscplx = {kk: np.iscomplexobj(lst0[kk]) for kk in dom.keys()} - - N = sum(iscplx.values()) * 2 + sum(not aa for aa in iscplx.values()) - n1 = int(np.ceil(np.sqrt(N))) - n2 = int(np.ceil(N / n1)) - - if n1 == n2 == 1: - fig = plt.figure(figsize=(UNIT / 2 * n1, UNIT / 2 * n2)) - axs = [plt.gca()] - else: - fig, axs = plt.subplots(n2, n1, figsize=(UNIT / 2 * n1, UNIT * 3 / 4 * n2)) - axs = list(axs.ravel()) - - for kk, cplx in iscplx.items(): - if isinstance(dom, ift.MultiDomain): - lstA = np.array([ss[kk].val for ss in lst]) - else: - lstA = np.array([ss.val for ss in lst]) - if cplx: - lstC = {"real": lstA.real, "imag": lstA.imag} - else: - lstC = {"": lstA} - for ll, lstD in lstC.items(): - axx = axs.pop(0) - lstD = lstD.ravel() - if lstD.size == 1: - axx.axvline(lstD[0]) - else: - axx.hist( - lstD.clip(-lim, lim - width), - alpha=0.5, - density=True, - bins=np.arange(-lim, lim + width, width), - ) - xs = np.linspace(-5, 5) - axx.set_yscale("log") - axx.plot(xs, np.exp(-(xs ** 2) / 2) / np.sqrt(2 * np.pi)) - redchisq = ( - np.sum((np.abs(lstD) if np.iscomplexobj(lstD) else lstD) ** 2) - / lstD.size - ) - axx.set_title(f"{kk} {ll} {redchisq:.1f}") - plt.tight_layout() - fig.savefig(fname) - plt.close(fig) - - -def _plot_spectra(state, op, name, directions): - fld = op.force(state.mean) - mi, ma = np.min(fld.val), np.max(fld.val) - - fdom, dom = fld.domain - freqs = np.array(fdom.coordinates) * 1e-6 - - plt.figure() - for indx, indy in directions: - plt.plot(freqs, fld.val[:, indx, indy]) - plt.ylim([mi, ma]) - plt.xlabel("MHz") - plt.tight_layout() - plt.savefig(name) - plt.close() - - -def _plot_mf(state, op, name, directions, movie_length): - # FIXME Write proper plotting routine which includes also directions and uncertainties - fld = op.force(state.mean) - - fdom, dom = fld.domain - freqs = np.array(fdom.coordinates) - nfreq = fld.shape[0] - my_asserteq(nfreq, len(freqs)) - mi, ma = np.min(fld.val), np.max(fld.val) - - name = splitext(name)[0] - for ii in range(nfreq): - ift.single_plot( - ift.makeField(dom, fld.val[ii]), - title=f"{freqs[ii]*1e-6:.0f} MHz", - vmin=mi, - vmax=ma, - cmap="inferno", - name=f"{name}_{ii:04.0f}.png", - ) - - if movie_length > 0: - framerate = nfreq / float(movie_length) - os.system( - f"ffmpeg -framerate {framerate} -i {name}_%04d.png -c:v libx264 -pix_fmt yuv420p -crf 23 -y {name}.mp4" - ) - for ii in range(nfreq): - os.remove(f"{name}_{ii:04.0f}.png") diff --git a/resolve/sky_model.py b/resolve/sky_model.py new file mode 100644 index 00000000..cf49309f --- /dev/null +++ b/resolve/sky_model.py @@ -0,0 +1,67 @@ +# SPDX-License-Identifier: GPL-3.0-or-later +# Copyright(C) 2021 Max-Planck-Society +# Author: Philipp Arras + +import nifty8 as ift + +from .constants import str2rad +from .points import PointInserter + + +def single_frequency_sky(cfg_section, list_point_sources=[]): + nx = cfg_section.getint("space npix x") + ny = cfg_section.getint("space npix y") + dx = str2rad(cfg_section["space fov x"]) / nx + dy = str2rad(cfg_section["space fov y"]) / ny + dom = ift.RGSpace([nx, ny], [dx, dy]) + + logdiffuse, cfm = cfm_from_cfg(cfg_section, {"space": dom}, "sky diffuse" ) + sky = logdiffuse.exp() + add = {"logdiffuse": logdiffuse, "logdiffuse power spectrum": cfm.power_spectrum} + + # Point sources + mode = cfg_section["point sources mode"] + if mode == "single": + ppos = [] + s = cfg_section["point sources locations"] + for xy in s.split(";"): + x, y = xy.split(",") + ppos.append((str2rad(x), str2rad(y))) + inserter = PointInserter(dom, ppos) + alpha = cfg_section.getfloat("point sources alpha") + q = cfg_section.getfloat("point sources q") + points = ift.InverseGammaOperator(inserter.domain, alpha=alpha, q=q/dom.scalar_dvol) + points = inserter @ points.ducktape("points") + sky = sky + points + add["points"] = points + + elif mode == "disable": + add["points"] = None + else: + raise ValueError(f"In order to disable point source component, set `point sources mode` to `disable`. Got: {mode}") + return sky, add + + +def cfm_from_cfg(cfg_section, domain_dct, prefix): + if not isinstance(domain_dct, dict): + domain_dct = {"": domain_dct} + cfm = ift.CorrelatedFieldMaker(_append_to_nonempty_string(prefix, " ")) + for key_prefix, dom in domain_dct.items(): + foo = _append_to_nonempty_string(key_prefix, " ") + kwargs = {kk: tuple(cfg_section.getfloat(f"{foo}{kk} {stat}") for stat in ["mean", "stddev"]) for kk in ["fluctuations", "loglogavgslope", "flexibility", "asperity"]} + kwargs["prefix"] = key_prefix + cfm.add_fluctuations(dom, **kwargs) + kwargs = { + "offset_mean": cfg_section.getfloat("zero mode offset"), + "offset_std": (cfg_section.getfloat("zero mode mean"), + cfg_section.getfloat("zero mode stddev")) + } + cfm.set_amplitude_total_offset(**kwargs) + op = cfm.finalize(prior_info=0) + return op, cfm + + +def _append_to_nonempty_string(s, append): + if s == "": + return s + return s + append -- GitLab From d05897c3174fcf09a5394ccec8d02b8ebc344215 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Tue, 9 Nov 2021 16:59:45 +0100 Subject: [PATCH 55/89] Try to fix MPI problem --- demo/imaging_with_automatic_weighting.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/demo/imaging_with_automatic_weighting.py b/demo/imaging_with_automatic_weighting.py index cf0f3038..5ce668d9 100644 --- a/demo/imaging_with_automatic_weighting.py +++ b/demo/imaging_with_automatic_weighting.py @@ -58,7 +58,7 @@ def main(cfg_file_name): assert obs.nfreq == obs.npol == 1 subcfg = cfg["weighting"] if subcfg["model"] == "cfm": - npix = 2500 + npix = 2500 # FIXME Use numbers from config file xs = np.log(obs.effective_uvwlen().val) xs -= np.min(xs) if not xs.shape[0] == xs.shape[2] == 1: @@ -112,6 +112,7 @@ def main(cfg_file_name): position = ift.MultiField.union([0.1*ift.from_random(sky.domain), position.extract(weightop.domain)]) if rve.mpi.mpi: + rve.mpi.comm.Barrier() if not master: position = None position = rve.mpi.comm.bcast(position, root=0) -- GitLab From cbadab38489aff9c68ca1e87c9fce3ab213bd61a Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Wed, 10 Nov 2021 12:12:30 +0100 Subject: [PATCH 56/89] Make wiener process more flexible --- resolve/multi_frequency/operators.py | 31 ++++++++++++++++++++-------- 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/resolve/multi_frequency/operators.py b/resolve/multi_frequency/operators.py index b77f98b2..20200fce 100644 --- a/resolve/multi_frequency/operators.py +++ b/resolve/multi_frequency/operators.py @@ -55,17 +55,30 @@ class WienerIntegrations(ift.LinearOperator): return ift.makeField(self._tgt(mode), res) -def IntWProcessInitialConditions(a0, b0, wfop): - for op in [a0, b0, wfop]: +def IntWProcessInitialConditions(a0, b0, wpop, irg_space=None): + for op in [a0, b0]: my_assert(ift.is_operator(op)) - my_asserteq(a0.target, b0.target, ift.makeDomain(wfop.target[1])) - bc = _FancyBroadcast(wfop.target) - factors = ift.full(wfop.target[0], 0) - factors = np.empty(wfop.target[0].shape) + + if ift.is_operator(wpop): + tgt = wpop.target + else: + tgt = irg_space, a0.target[0] + my_asserteq(a0.target[0], b0.target[0], tgt[1]) + + sdom = tgt[0] + + bc = _FancyBroadcast(tgt) + factors = ift.full(sdom, 0) + factors = np.empty(sdom.shape) factors[0] = 0 - factors[1:] = np.cumsum(wfop.target[0].dvol) - factors = ift.makeField(wfop.target[0], factors) - return wfop + bc @ a0 + ift.DiagonalOperator(factors, wfop.target, 0) @ bc @ b0 + factors[1:] = np.cumsum(sdom.dvol) + factors = ift.makeField(sdom, factors) + res = bc @ a0 + ift.DiagonalOperator(factors, tgt, 0) @ bc @ b0 + + if wpop is None: + return res + else: + return wpop + res class _FancyBroadcast(ift.LinearOperator): -- GitLab From 2ba87495d6d276d738dcb776d19355824dc75a1e Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Wed, 10 Nov 2021 15:14:05 +0100 Subject: [PATCH 57/89] Make Wiener integration more flexible --- resolve/multi_frequency/operators.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/resolve/multi_frequency/operators.py b/resolve/multi_frequency/operators.py index 20200fce..af0ef524 100644 --- a/resolve/multi_frequency/operators.py +++ b/resolve/multi_frequency/operators.py @@ -27,7 +27,9 @@ class WienerIntegrations(ift.LinearOperator): self._target = ift.makeDomain((time_domain, remaining_domain)) dom = ift.UnstructuredDomain((2, time_domain.size - 1)), remaining_domain self._domain = ift.makeDomain(dom) - self._volumes = time_domain.dvol[:, None, None] + self._volumes = time_domain.dvol + for _ in range(len(remaining_domain.shape)): + self._volumes = self._volumes[..., np.newaxis] self._capability = self.TIMES | self.ADJOINT_TIMES def apply(self, x, mode): -- GitLab From ae59b63770c09f9cd647b9217eef5c9a16973617 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Mon, 15 Nov 2021 17:30:48 +0100 Subject: [PATCH 58/89] Work on standard reconstruction script --- demo/imaging_with_automatic_weighting.py | 108 ++++++++++++----------- 1 file changed, 57 insertions(+), 51 deletions(-) diff --git a/demo/imaging_with_automatic_weighting.py b/demo/imaging_with_automatic_weighting.py index 5ce668d9..8ca87cec 100644 --- a/demo/imaging_with_automatic_weighting.py +++ b/demo/imaging_with_automatic_weighting.py @@ -80,61 +80,66 @@ def main(cfg_file_name): full_lh = rve.ImagingLikelihood(obs, sky, inverse_covariance_operator=weightop) position = 0.1 * ift.from_random(full_lh.domain) - common = {"output_directory": cfg["output"]["directory"], - "plottable_operators": operators, - "overwrite": True} - if master: - common["comm"] = rve.mpi.comm_self - if enable_points: - # Points only, MAP - lh = rve.ImagingLikelihood(obs, operators["points"]) - mini = ift.NewtonCG(ift.GradientNormController(name="hamiltonian", iteration_limit=4)) - sl = ift.optimize_kl(lh, 1, 0, mini, None, None, initial_position=position, **common) - position = sl.local_item(0) - - # Points constant, diffuse only, MAP - lh = rve.ImagingLikelihood(obs, sky) - mini = ift.NewtonCG(ift.GradientNormController(name="hamiltonian", iteration_limit=20)) - sl = ift.optimize_kl(lh, 1, 0, mini, None, None, initial_position=position, initial_index=1, **common) - position = sl.local_item(0) - - # Only weighting, MGVI - cst = sky.domain.keys() - mini = ift.VL_BFGS(ift.GradientNormController(name="bfgs", iteration_limit=20)) - ic_sampling = ift.AbsDeltaEnergyController(0.1, 3, 100, name="Sampling") - sl = ift.optimize_kl(full_lh, 5, 3, mini, ic_sampling, None, - constants=cst, point_estimates=cst, - initial_position=position, initial_index=2, **common) - position = sl.local_item(0) - - # Reset sky - position = ift.MultiField.union([0.1*ift.from_random(sky.domain), position.extract(weightop.domain)]) - - if rve.mpi.mpi: - rve.mpi.comm.Barrier() - if not master: - position = None - position = rve.mpi.comm.bcast(position, root=0) - ift.random.push_sseq_from_seed(42) - common["comm"] = rve.mpi.comm + common = {"plottable_operators": operators, "overwrite": True} + def get_mini(iglobal): - if iglobal < 7: - return ift.NewtonCG(ift.GradientNormController(name="kl", iteration_limit=15)) - elif iglobal < 7 + 5: + if iglobal == 0: + return ift.NewtonCG(ift.GradientNormController(name="hamiltonian", iteration_limit=4)) + if iglobal == 1: + return ift.NewtonCG(ift.GradientNormController(name="hamiltonian", iteration_limit=20)) + return ift.VL_BFGS(ift.GradientNormController(name="bfgs", iteration_limit=20)) + + def get_sampling(iglobal): + if iglobal in [0, 1]: + return None + return ift.AbsDeltaEnergyController(deltaE=0.5, convergence_level=3, iteration_limit=500, + name="Sampling") + + def get_cst(iglobal): + if iglobal in [0, 1]: + return [] + return sky.domain.keys() + + def get_lh(iglobal): + if iglobal == 0: + return rve.ImagingLikelihood(obs, operators["points"]) + if iglobal == 1: + return rve.ImagingLikelihood(obs, sky) + return full_lh + + def get_n_samples(iglobal): + if iglobal in [0, 1]: + return 0 + return 3 + + def get_comm(iglobal): + return rve.mpi.comm + + def callback(sl, iglobal): + lh = get_lh(iglobal) + s = ift.extra.minisanity(lh.data, lh.metric_at_pos, lh.model_data, sl) + if rve.mpi.master: + print(s) + + sl = ift.optimize_kl(get_lh, 7, get_n_samples, get_mini, get_sampling, None, + constants=get_cst, point_estimates=get_cst, initial_position=position, + comm=get_comm, callback=callback, + output_directory=cfg["output"]["directory"] + "initial", **common) + + # Reset sky + position = ift.MultiField.union([0.1*ift.from_random(sky.domain), + position.extract(weightop.domain)]) + + def get_mini(iglobal): + if iglobal < 5: return ift.VL_BFGS(ift.GradientNormController(name="kl", iteration_limit=15)) - else: - return ift.NewtonCG(ift.GradientNormController(name="kl", iteration_limit=15)) + return ift.NewtonCG(ift.GradientNormController(name="kl", iteration_limit=15)) def get_sampling(iglobal): - if iglobal < 7: - lim = 200 - elif iglobal < 7 + 5: - lim = 700 - else: - lim = 700 - return ift.AbsDeltaEnergyController(deltaE=0.5, iteration_limit=lim, name="Sampling") + return ift.AbsDeltaEnergyController(deltaE=0.5, convergence_level=3, iteration_limit=100, + name="Sampling") def get_cst(iglobal): res = [] @@ -144,10 +149,11 @@ def main(cfg_file_name): res += list(operators["points"].domain.keys()) return res - # First MGVI diffuse sky, then sky + weighting simultaneously sl = ift.optimize_kl(full_lh, 35, 6, get_mini, get_sampling, None, constants=get_cst, point_estimates=get_cst, - initial_position=position, initial_index=7, **common) + initial_position=position, comm=rve.mpi.comm, + output_directory=cfg["output"]["directory"], **common) + if __name__ == "__main__": -- GitLab From 7af07bb77bb98d2321d3a7786d86a3bd8655de36 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Mon, 15 Nov 2021 19:40:05 +0100 Subject: [PATCH 59/89] Use ducc version from PyPI --- Dockerfile | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Dockerfile b/Dockerfile index 36aca4e1..66bfad0f 100644 --- a/Dockerfile +++ b/Dockerfile @@ -6,8 +6,7 @@ RUN apt-get update -qq && apt-get install -qq git # Actual dependencies RUN apt-get update -qq && apt-get install -qq python3-pip casacore-dev python3-matplotlib -RUN pip3 install scipy git+https://gitlab.mpcdf.mpg.de/ift/nifty.git@NIFTy_8 -RUN pip3 install git+https://gitlab.mpcdf.mpg.de/mtr/ducc.git@ducc0 +RUN pip3 install scipy git+https://gitlab.mpcdf.mpg.de/ift/nifty.git@NIFTy_8 ducc0 # Optional dependencies RUN pip3 install astropy RUN apt-get install -qq python3-mpi4py -- GitLab From 420f17622dd79eb4d4029d3d6786ed5d1adc89f5 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Mon, 15 Nov 2021 19:40:18 +0100 Subject: [PATCH 60/89] Cosmetics --- demo/imaging_with_automatic_weighting.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/demo/imaging_with_automatic_weighting.py b/demo/imaging_with_automatic_weighting.py index 8ca87cec..d91c8c9a 100644 --- a/demo/imaging_with_automatic_weighting.py +++ b/demo/imaging_with_automatic_weighting.py @@ -83,7 +83,6 @@ def main(cfg_file_name): common = {"plottable_operators": operators, "overwrite": True} - def get_mini(iglobal): if iglobal == 0: return ift.NewtonCG(ift.GradientNormController(name="hamiltonian", iteration_limit=4)) @@ -126,7 +125,7 @@ def main(cfg_file_name): sl = ift.optimize_kl(get_lh, 7, get_n_samples, get_mini, get_sampling, None, constants=get_cst, point_estimates=get_cst, initial_position=position, comm=get_comm, callback=callback, - output_directory=cfg["output"]["directory"] + "initial", **common) + output_directory=cfg["output"]["directory"] + "_initial", **common) # Reset sky position = ift.MultiField.union([0.1*ift.from_random(sky.domain), @@ -153,7 +152,6 @@ def main(cfg_file_name): constants=get_cst, point_estimates=get_cst, initial_position=position, comm=rve.mpi.comm, output_directory=cfg["output"]["directory"], **common) - if __name__ == "__main__": -- GitLab From b2e7bc9a1a92521718373c2485ddffcc51a0d758 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Tue, 16 Nov 2021 16:55:13 +0100 Subject: [PATCH 61/89] IRGSpace: symmetrical pixels --- resolve/multi_frequency/irg_space.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/resolve/multi_frequency/irg_space.py b/resolve/multi_frequency/irg_space.py index a2a23bc7..d7a1f7de 100644 --- a/resolve/multi_frequency/irg_space.py +++ b/resolve/multi_frequency/irg_space.py @@ -59,10 +59,13 @@ class IRGSpace(ift.StructuredDomain): @property def dvol(self): - """This has one pixel less than :meth:`shape`. Not sure if this is - okay. - """ - return np.diff(np.array(self._coordinates)) + """Assume that the coordinates are the center of symmetric pixels.""" + c = np.array(self._coordinates) + bounds = np.empty(self.size + 1) + bounds[1:-1] = c[:-1] + 0.5*np.diff(c) + bounds[0] = c[0] - 0.5*(c[1] - c[0]) + bounds[-1] = c[-1] + 0.5*(c[-1] - c[-2]) + return np.diff(bounds) @property def coordinates(self): -- GitLab From b33d94709affcce1a13be4b88b5075f8de464aab Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 25 Nov 2021 15:38:00 +0100 Subject: [PATCH 62/89] IRGSpace: Introduce distances --- resolve/multi_frequency/irg_space.py | 4 ++++ resolve/multi_frequency/operators.py | 8 +++----- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/resolve/multi_frequency/irg_space.py b/resolve/multi_frequency/irg_space.py index d7a1f7de..34ac470c 100644 --- a/resolve/multi_frequency/irg_space.py +++ b/resolve/multi_frequency/irg_space.py @@ -67,6 +67,10 @@ class IRGSpace(ift.StructuredDomain): bounds[-1] = c[-1] + 0.5*(c[-1] - c[-2]) return np.diff(bounds) + @property + def distances(self): + return np.diff(self._coordinates) + @property def coordinates(self): return self._coordinates diff --git a/resolve/multi_frequency/operators.py b/resolve/multi_frequency/operators.py index af0ef524..a84940cb 100644 --- a/resolve/multi_frequency/operators.py +++ b/resolve/multi_frequency/operators.py @@ -27,7 +27,7 @@ class WienerIntegrations(ift.LinearOperator): self._target = ift.makeDomain((time_domain, remaining_domain)) dom = ift.UnstructuredDomain((2, time_domain.size - 1)), remaining_domain self._domain = ift.makeDomain(dom) - self._volumes = time_domain.dvol + self._volumes = time_domain.distances for _ in range(len(remaining_domain.shape)): self._volumes = self._volumes[..., np.newaxis] self._capability = self.TIMES | self.ADJOINT_TIMES @@ -42,9 +42,7 @@ class WienerIntegrations(ift.LinearOperator): x = x.val res = np.zeros(self._target.shape) res[from_second] = np.cumsum(x[second], axis=0) - res[from_second] = ( - res[from_second] + res[no_border] - ) / 2 * self._volumes + x[first] + res[from_second] = (res[from_second] + res[no_border]) / 2 * self._volumes + x[first] res[from_second] = np.cumsum(res[from_second], axis=0) else: x = x.val_rw() @@ -73,7 +71,7 @@ def IntWProcessInitialConditions(a0, b0, wpop, irg_space=None): factors = ift.full(sdom, 0) factors = np.empty(sdom.shape) factors[0] = 0 - factors[1:] = np.cumsum(sdom.dvol) + factors[1:] = np.cumsum(sdom.distances) factors = ift.makeField(sdom, factors) res = bc @ a0 + ift.DiagonalOperator(factors, tgt, 0) @ bc @ b0 -- GitLab From eac2293a2e31b1061eb2299601f15e7e23b78139 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 25 Nov 2021 15:38:16 +0100 Subject: [PATCH 63/89] Cleanup --- resolve/__init__.py | 1 - resolve/multi_frequency/operators.py | 30 ---------------------------- test/test_general.py | 11 ---------- 3 files changed, 42 deletions(-) diff --git a/resolve/__init__.py b/resolve/__init__.py index 8321f28e..194bd883 100644 --- a/resolve/__init__.py +++ b/resolve/__init__.py @@ -15,7 +15,6 @@ from .data.ms_import import * from .multi_frequency.irg_space import IRGSpace from .multi_frequency.operators import ( IntWProcessInitialConditions, - MfWeightingInterpolation, WienerIntegrations, ) from .data.observation import * diff --git a/resolve/multi_frequency/operators.py b/resolve/multi_frequency/operators.py index a84940cb..8d415067 100644 --- a/resolve/multi_frequency/operators.py +++ b/resolve/multi_frequency/operators.py @@ -5,7 +5,6 @@ import nifty8 as ift import numpy as np -from ..polarization_space import PolarizationSpace from ..util import my_assert, my_assert_isinstance, my_asserteq from .irg_space import IRGSpace @@ -96,32 +95,3 @@ class _FancyBroadcast(ift.LinearOperator): else: res = np.sum(x.val, axis=0) return ift.makeField(self._tgt(mode), res) - - -class MfWeightingInterpolation(ift.LinearOperator): - def __init__(self, eff_uvw, domain): - domain = ift.DomainTuple.make(domain) - my_asserteq(domain.shape[0], eff_uvw.shape[2]) # freqaxis - self._domain = domain - nrow, nfreq = eff_uvw.shape[1:] - tgt = [PolarizationSpace("I")] + [ift.UnstructuredDomain(aa) for aa in [nrow, nfreq]] - self._target = ift.DomainTuple.make(tgt) - self._capability = self.TIMES | self.ADJOINT_TIMES - # FIXME Try to unify all those operators which loop over freq dimension - self._ops = [] - for ii in range(nfreq): - op = ift.LinearInterpolator(domain[1], eff_uvw.val[:, :, ii]) - self._ops.append(op) - my_asserteq(self.target.shape[0], 1) - - def apply(self, x, mode): - self._check_input(x, mode) - res = np.empty(self._tgt(mode).shape) - if mode == self.TIMES: - for ii, op in enumerate(self._ops): - res[0, :, ii] = op(ift.makeField(op.domain, x.val[ii])).val - else: - for ii, op in enumerate(self._ops): - op = op.adjoint - res[ii] = op(ift.makeField(op.domain, x.val[0, :, ii])).val - return ift.makeField(self._tgt(mode), res) diff --git a/test/test_general.py b/test/test_general.py index 4936fb2f..7cc52b79 100644 --- a/test/test_general.py +++ b/test/test_general.py @@ -311,17 +311,6 @@ def test_randomonmaster(): finvalid() -def test_mfweighting(): - nrow = 100 - nchan = 4 - npol = 1 - effuv = ift.random.current_rng().random((npol, nrow, nchan)) - effuv = ift.makeField([ift.UnstructuredDomain(ii) for ii in effuv.shape], effuv) - dom = ift.UnstructuredDomain(nchan), ift.RGSpace(npix, 2 * np.max(effuv.val) / npix) - op = rve.MfWeightingInterpolation(effuv, dom) - ift.extra.check_linear_operator(op) - - def test_mf_response(): ms = join(direc, "CYG-D-6680-64CH-10S.ms") obs = rve.ms2observations(ms, "DATA", False, 0, "stokesiavg")[0] -- GitLab From c9ca50c0b2799bea32d85b73c770b60fc4b47c27 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Wed, 17 Nov 2021 11:32:41 +0100 Subject: [PATCH 64/89] Add convenience functions for calibration --- resolve/calibration.py | 2 +- resolve/data/observation.py | 7 +++++++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/resolve/calibration.py b/resolve/calibration.py index bcd322c8..a9475581 100644 --- a/resolve/calibration.py +++ b/resolve/calibration.py @@ -34,7 +34,7 @@ class CalibrationDistributor(ift.LinearOperator): my_assert(np.issubdtype(ant_col.dtype, np.integer)) my_assert(set(np.unique(ant_col)) <= set(antenna_dct.keys())) my_assert(np.min(ant_col) >= 0) - my_assert(np.max(ant_col) < len(antenna_dct)) + my_assert(np.max(ant_col) <= max(antenna_dct.keys())) self._domain = ift.DomainTuple.make(domain) self._target = ift.DomainTuple.make(target) my_asserteq(len(self._domain), 4) diff --git a/resolve/data/observation.py b/resolve/data/observation.py index d3f13ea6..3c3ec0f0 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -597,6 +597,13 @@ class Observation(BaseObservation): print(f"INFO: Flag baseline {ant1_name}-{ant2_name}, {np.sum(ind)}/{self.nrow} rows flagged.") return Observation(self._antpos, self._vis, wgt, self._polarization, self._freq, self._auxiliary_tables) + @property + def nbaselines(self): + return len(self.baselines()) + + def baselines(self): + return set((a1, a2) for a1, a2 in zip(self.ant1, self.ant2)) + @property def uvw(self): return self._antpos.uvw -- GitLab From 14979a1096772698e771119e616be9922e53ded3 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Wed, 17 Nov 2021 15:15:06 +0100 Subject: [PATCH 65/89] IRGSpace: add binbounds() --- resolve/multi_frequency/irg_space.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/resolve/multi_frequency/irg_space.py b/resolve/multi_frequency/irg_space.py index 34ac470c..e67ee8d1 100644 --- a/resolve/multi_frequency/irg_space.py +++ b/resolve/multi_frequency/irg_space.py @@ -60,12 +60,15 @@ class IRGSpace(ift.StructuredDomain): @property def dvol(self): """Assume that the coordinates are the center of symmetric pixels.""" + return np.diff(self.binbounds()) + + def binbounds(self): c = np.array(self._coordinates) bounds = np.empty(self.size + 1) bounds[1:-1] = c[:-1] + 0.5*np.diff(c) bounds[0] = c[0] - 0.5*(c[1] - c[0]) bounds[-1] = c[-1] + 0.5*(c[-1] - c[-2]) - return np.diff(bounds) + return bounds @property def distances(self): -- GitLab From bfb48146b75a21e9b004b5fb05f5169358f04902 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Tue, 23 Nov 2021 11:48:24 +0100 Subject: [PATCH 66/89] Cygnus A demo: finalize --- demo/cygnusa.cfg | 3 +- demo/imaging_with_automatic_weighting.py | 37 ++++++++++++++---------- 2 files changed, 24 insertions(+), 16 deletions(-) diff --git a/demo/cygnusa.cfg b/demo/cygnusa.cfg index c9f8bb57..6d7e0e1c 100644 --- a/demo/cygnusa.cfg +++ b/demo/cygnusa.cfg @@ -36,7 +36,8 @@ point sources q = 0.2 enable = True model = cfm -npix = 2500 +npix = 1000 +zeropadding factor = 2 zero mode offset = 0 zero mode mean = 2 zero mode stddev = 2 diff --git a/demo/imaging_with_automatic_weighting.py b/demo/imaging_with_automatic_weighting.py index d91c8c9a..5fcf278f 100644 --- a/demo/imaging_with_automatic_weighting.py +++ b/demo/imaging_with_automatic_weighting.py @@ -58,13 +58,17 @@ def main(cfg_file_name): assert obs.nfreq == obs.npol == 1 subcfg = cfg["weighting"] if subcfg["model"] == "cfm": - npix = 2500 # FIXME Use numbers from config file + import ducc0 + + npix = subcfg.getint("npix") + fac = subcfg.getfloat("zeropadding factor") + npix_padded = ducc0.fft.good_size(int(np.round(npix*fac))) + xs = np.log(obs.effective_uvwlen().val) xs -= np.min(xs) if not xs.shape[0] == xs.shape[2] == 1: raise RuntimeError - # FIXME Use Integrated Wiener process - dom = ift.RGSpace(npix, 2 * np.max(xs) / npix) + dom = ift.RGSpace(npix_padded, np.max(xs) / npix) logwgt, cfm = rve.cfm_from_cfg(subcfg, dom, "invcov") li = ift.LinearInterpolator(dom, xs[0].T) conv = rve.DomainChangerAndReshaper(li.target, obs.weight.domain) @@ -81,7 +85,7 @@ def main(cfg_file_name): full_lh = rve.ImagingLikelihood(obs, sky, inverse_covariance_operator=weightop) position = 0.1 * ift.from_random(full_lh.domain) - common = {"plottable_operators": operators, "overwrite": True} + common = {"plottable_operators": operators, "overwrite": True, "return_final_position": True} def get_mini(iglobal): if iglobal == 0: @@ -122,14 +126,17 @@ def main(cfg_file_name): if rve.mpi.master: print(s) - sl = ift.optimize_kl(get_lh, 7, get_n_samples, get_mini, get_sampling, None, - constants=get_cst, point_estimates=get_cst, initial_position=position, - comm=get_comm, callback=callback, - output_directory=cfg["output"]["directory"] + "_initial", **common) + if True: + _, position = ift.optimize_kl(get_lh, 7, get_n_samples, get_mini, get_sampling, None, + constants=get_cst, point_estimates=get_cst, + initial_position=position, comm=get_comm, callback=callback, + output_directory=cfg["output"]["directory"] + "_initial", + **common) + else: + position = ift.ResidualSampleList.load_mean("Cygnus_A_2052MHz_initial/pickle/last") - # Reset sky - position = ift.MultiField.union([0.1*ift.from_random(sky.domain), - position.extract(weightop.domain)]) + # Reset diffuse component + position = ift.MultiField.union([position, 0.1*ift.from_random(operators["logdiffuse"].domain)]) def get_mini(iglobal): if iglobal < 5: @@ -148,10 +155,10 @@ def main(cfg_file_name): res += list(operators["points"].domain.keys()) return res - sl = ift.optimize_kl(full_lh, 35, 6, get_mini, get_sampling, None, - constants=get_cst, point_estimates=get_cst, - initial_position=position, comm=rve.mpi.comm, - output_directory=cfg["output"]["directory"], **common) + ift.optimize_kl(full_lh, 35, 6, get_mini, get_sampling, None, + constants=get_cst, point_estimates=get_cst, + initial_position=position, comm=rve.mpi.comm, + output_directory=cfg["output"]["directory"], **common) if __name__ == "__main__": -- GitLab From 98c4fd4adaa4e0d9c5f29de0fa86fa10f7c97ed7 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 26 Nov 2021 11:24:03 +0100 Subject: [PATCH 67/89] Observation: Do not enforce single-precision anymore --- resolve/data/observation.py | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/resolve/data/observation.py b/resolve/data/observation.py index 3c3ec0f0..2da44d84 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -284,14 +284,8 @@ class Observation(BaseObservation): my_asserteq(weight.shape, vis.shape) my_asserteq(vis.shape, (len(polarization), nrows, len(freq))) my_asserteq(nrows, vis.shape[1]) - if vis.dtype != np.complex64: - print(f"Warning: vis.dtype is {vis.dtype}. Casting to np.complex64") - vis = vis.astype(np.complex64) - if weight.dtype != np.float32: - print(f"Warning: weight.dtype is {weight.dtype}. Casting to np.float32") - weight = weight.astype(np.float32) - my_asserteq(vis.dtype, np.complex64) - my_asserteq(weight.dtype, np.float32) + #my_asserteq(vis.dtype, np.complex64) + #my_asserteq(weight.dtype, np.float32) my_assert(np.all(weight >= 0.0)) my_assert(np.all(np.isfinite(vis[weight > 0.]))) my_assert(np.all(np.isfinite(weight))) -- GitLab From 3afd8f125c81636d90eef02e94508dbb0cf8d58c Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 26 Nov 2021 17:01:32 +0100 Subject: [PATCH 68/89] IRGSpace: binbounds for single bin (-inf, inf) --- resolve/multi_frequency/irg_space.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/resolve/multi_frequency/irg_space.py b/resolve/multi_frequency/irg_space.py index e67ee8d1..a6c3d4bd 100644 --- a/resolve/multi_frequency/irg_space.py +++ b/resolve/multi_frequency/irg_space.py @@ -63,6 +63,8 @@ class IRGSpace(ift.StructuredDomain): return np.diff(self.binbounds()) def binbounds(self): + if len(self._coordinates) == 1: + return np.array([-np.inf, np.inf]) c = np.array(self._coordinates) bounds = np.empty(self.size + 1) bounds[1:-1] = c[:-1] + 0.5*np.diff(c) -- GitLab From 94242fdfb3dc091c11fcc260da2411c05240603a Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Mon, 29 Nov 2021 16:40:45 +0100 Subject: [PATCH 69/89] IRGSpace: Add coordinates to __repr__ --- resolve/multi_frequency/irg_space.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/resolve/multi_frequency/irg_space.py b/resolve/multi_frequency/irg_space.py index a6c3d4bd..bee19fdf 100644 --- a/resolve/multi_frequency/irg_space.py +++ b/resolve/multi_frequency/irg_space.py @@ -38,7 +38,7 @@ class IRGSpace(ift.StructuredDomain): self._coordinates = tuple(bb) def __repr__(self): - return f"IRGSpace(shape={self.shape}, coordinates=...)" + return f"IRGSpace(shape={self.shape}, coordinates={self._coordinates})" @property def harmonic(self): -- GitLab From a7cfac4cfaef7c20bdb599344f9cbf21a2c078ab Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 25 Nov 2021 15:41:36 +0100 Subject: [PATCH 70/89] Restructure --- resolve/__init__.py | 4 ++-- .../operators.py => integrated_wiener_process.py} | 0 resolve/{multi_frequency => }/irg_space.py | 0 resolve/multi_frequency/__init__.py | 0 resolve/response.py | 2 +- 5 files changed, 3 insertions(+), 3 deletions(-) rename resolve/{multi_frequency/operators.py => integrated_wiener_process.py} (100%) rename resolve/{multi_frequency => }/irg_space.py (100%) delete mode 100644 resolve/multi_frequency/__init__.py diff --git a/resolve/__init__.py b/resolve/__init__.py index 194bd883..9613ce36 100644 --- a/resolve/__init__.py +++ b/resolve/__init__.py @@ -12,8 +12,8 @@ from .mosaicing import * from .mpi import onlymaster from .mpi_operators import * from .data.ms_import import * -from .multi_frequency.irg_space import IRGSpace -from .multi_frequency.operators import ( +from .irg_space import IRGSpace +from .integrated_wiener_process import ( IntWProcessInitialConditions, WienerIntegrations, ) diff --git a/resolve/multi_frequency/operators.py b/resolve/integrated_wiener_process.py similarity index 100% rename from resolve/multi_frequency/operators.py rename to resolve/integrated_wiener_process.py diff --git a/resolve/multi_frequency/irg_space.py b/resolve/irg_space.py similarity index 100% rename from resolve/multi_frequency/irg_space.py rename to resolve/irg_space.py diff --git a/resolve/multi_frequency/__init__.py b/resolve/multi_frequency/__init__.py deleted file mode 100644 index e69de29b..00000000 diff --git a/resolve/response.py b/resolve/response.py index c9271339..14e0caff 100644 --- a/resolve/response.py +++ b/resolve/response.py @@ -13,7 +13,7 @@ import nifty8 as ift from .constants import SPEEDOFLIGHT from .global_config import epsilon, nthreads, wgridding, verbosity, np_dtype -from .multi_frequency.irg_space import IRGSpace +from .irg_space import IRGSpace from .data.observation import Observation from .util import my_assert, my_assert_isinstance, my_asserteq -- GitLab From 1d062cb3de0313ef57b21324e8f638786cd0665d Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 25 Nov 2021 15:42:57 +0100 Subject: [PATCH 71/89] Cleanup --- resolve/sky_model.py | 2 +- test/test_general.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/resolve/sky_model.py b/resolve/sky_model.py index cf49309f..82ae12c4 100644 --- a/resolve/sky_model.py +++ b/resolve/sky_model.py @@ -8,7 +8,7 @@ from .constants import str2rad from .points import PointInserter -def single_frequency_sky(cfg_section, list_point_sources=[]): +def single_frequency_sky(cfg_section): nx = cfg_section.getint("space npix x") ny = cfg_section.getint("space npix y") dx = str2rad(cfg_section["space fov x"]) / nx diff --git a/test/test_general.py b/test/test_general.py index 7cc52b79..3405b254 100644 --- a/test/test_general.py +++ b/test/test_general.py @@ -15,7 +15,6 @@ np.seterr(all="raise") direc = "/data/" nthreads = 1 -ift.set_nthreads(nthreads) rve.set_nthreads(nthreads) rve.set_epsilon(1e-4) rve.set_wgridding(False) -- GitLab From 8b09d07093a51ff1b4f1c017c782815b82b3d169 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 25 Nov 2021 15:53:08 +0100 Subject: [PATCH 72/89] fixup! Restructure --- resolve/integrated_wiener_process.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/resolve/integrated_wiener_process.py b/resolve/integrated_wiener_process.py index 8d415067..c60eb62a 100644 --- a/resolve/integrated_wiener_process.py +++ b/resolve/integrated_wiener_process.py @@ -5,8 +5,8 @@ import nifty8 as ift import numpy as np -from ..util import my_assert, my_assert_isinstance, my_asserteq from .irg_space import IRGSpace +from .util import my_assert, my_assert_isinstance, my_asserteq class WienerIntegrations(ift.LinearOperator): -- GitLab From 442a8c8d4ae03d0840810938f60f1ff850a98c03 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 25 Nov 2021 16:31:23 +0100 Subject: [PATCH 73/89] Start general interface for mf skys --- resolve/sky_model.py | 58 ++++++++++++++++++++++++++++++++++++----- resolve/util.py | 34 ++++++++++++++++++++++-- test/mf_sky.cfg | 37 ++++++++++++++++++++++++++ test/test_sky_models.py | 18 +++++++++++++ 4 files changed, 139 insertions(+), 8 deletions(-) create mode 100644 test/mf_sky.cfg create mode 100644 test/test_sky_models.py diff --git a/resolve/sky_model.py b/resolve/sky_model.py index 82ae12c4..5751a718 100644 --- a/resolve/sky_model.py +++ b/resolve/sky_model.py @@ -3,17 +3,18 @@ # Author: Philipp Arras import nifty8 as ift +import numpy as np from .constants import str2rad +from .irg_space import IRGSpace from .points import PointInserter +from .polarization_space import PolarizationSpace +from .simple_operators import DomainChangerAndReshaper +from .util import assert_sky_domain def single_frequency_sky(cfg_section): - nx = cfg_section.getint("space npix x") - ny = cfg_section.getint("space npix y") - dx = str2rad(cfg_section["space fov x"]) / nx - dy = str2rad(cfg_section["space fov y"]) / ny - dom = ift.RGSpace([nx, ny], [dx, dy]) + dom = _spatial_dom(cfg_section) logdiffuse, cfm = cfm_from_cfg(cfg_section, {"space": dom}, "sky diffuse" ) sky = logdiffuse.exp() @@ -42,13 +43,41 @@ def single_frequency_sky(cfg_section): return sky, add +def multi_frequency_sky(cfg): + if "point sources mode" in cfg: + raise NotImplementedError("Point sources are not supported yet.") + + sdom = _spatial_dom(cfg) + + if cfg["freq mode"] == "cfm": + fnpix, df = cfg.getfloat("freq npix"), cfg.getfloat("freq pixel size") + fdom = IRGSpace(cfg.getfloat("freq start") + np.arange(fnpix)*df) + fdom_rg = ift.RGSpace(fnpix, df) + logdiffuse, cfm = cfm_from_cfg(cfg, {"freq": fdom_rg, "space": sdom}, "sky diffuse") + sky = logdiffuse.exp() + add = {"logdiffuse": logdiffuse} + fampl, sampl = list(cfm.get_normalized_amplitudes()) + add["freq normalized power spectrum"] = fampl**2 + add["space normalized power spectrum"] = sampl**2 + + tgt = PolarizationSpace("I"), IRGSpace([np.nan]), fdom, sdom + sky = DomainChangerAndReshaper(sky.target, tgt) @ sky + elif cfg["freq mode"] == "integrated wiener process": + raise NotImplementedError + else: + raise RuntimeError + + assert_sky_domain(sky.target) + return sky, add + + def cfm_from_cfg(cfg_section, domain_dct, prefix): if not isinstance(domain_dct, dict): domain_dct = {"": domain_dct} cfm = ift.CorrelatedFieldMaker(_append_to_nonempty_string(prefix, " ")) for key_prefix, dom in domain_dct.items(): foo = _append_to_nonempty_string(key_prefix, " ") - kwargs = {kk: tuple(cfg_section.getfloat(f"{foo}{kk} {stat}") for stat in ["mean", "stddev"]) for kk in ["fluctuations", "loglogavgslope", "flexibility", "asperity"]} + kwargs = {kk: _parse_or_none(cfg_section, f"{foo}{kk}") for kk in ["fluctuations", "loglogavgslope", "flexibility", "asperity"]} kwargs["prefix"] = key_prefix cfm.add_fluctuations(dom, **kwargs) kwargs = { @@ -65,3 +94,20 @@ def _append_to_nonempty_string(s, append): if s == "": return s return s + append + + +def _spatial_dom(cfg): + nx = cfg.getint("space npix x") + ny = cfg.getint("space npix y") + dx = str2rad(cfg["space fov x"]) / nx + dy = str2rad(cfg["space fov y"]) / ny + return ift.RGSpace([nx, ny], [dx, dy]) + + +def _parse_or_none(cfg, key): + key0 = f"{key} mean" + key1 = f"{key} stddev" + if cfg[key0] == "None" and cfg[key1] == "None": + return None + if key0 in cfg: + return (cfg.getfloat(key0), cfg.getfloat(key1)) diff --git a/resolve/util.py b/resolve/util.py index 8093ef2e..9f015fdc 100644 --- a/resolve/util.py +++ b/resolve/util.py @@ -3,9 +3,8 @@ # Author: Philipp Arras import matplotlib.pyplot as plt -import numpy as np - import nifty8 as ift +import numpy as np def my_assert(*conds): @@ -124,3 +123,34 @@ class DomainChanger(ift.LinearOperator): def apply(self, x, mode): return ift.makeField(self._tgt(mode), x.val) + + +def assert_sky_domain(dom): + """Check that input fulfils resolve's conventions of a sky domain. + + A sky domain is a DomainTuple: + + dom = (pdom, tdom, fdom, sdom) + + where `pdom` is a `PolarizationSpace`, `tdom` and `fdom` are an `IRGSpace`, +and `sdom` is a two-dimensional `RGSpace`. + + Parameters + ---------- + dom : DomainTuple + Object that is checked to fulfil the properties. + """ + from .irg_space import IRGSpace + from .polarization_space import PolarizationSpace + + my_assert_isinstance(dom, ift.DomainTuple) + my_asserteq(len(dom), 4) + pdom, tdom, fdom, sdom = dom + for ii in range(3): + my_asserteq(len(dom[ii].shape), 1) + my_asserteq(len(sdom.shape), 2) + my_assert_isinstance(pdom, PolarizationSpace) + my_assert_isinstance(tdom, IRGSpace) + my_assert_isinstance(fdom, IRGSpace) + my_assert_isinstance(sdom, ift.RGSpace) + my_asserteq(len(sdom.shape), 2) diff --git a/test/mf_sky.cfg b/test/mf_sky.cfg new file mode 100644 index 00000000..d486869e --- /dev/null +++ b/test/mf_sky.cfg @@ -0,0 +1,37 @@ +[sky] +space npix x = 500 +space npix y = 500 +space fov x = 1000muas +space fov y = 1000muas + +# point sources mode = single +# point sources locations = 0deg,0deg;0.35as,-0.22as +# point sources alpha = 0.5 +# point sources q = 0.2 + +zero mode offset = 0 +zero mode mean = 1 +zero mode stddev = 0.3 + +space fluctuations mean = 1.5 +space fluctuations stddev = 1. +space loglogavgslope mean = -3 +space loglogavgslope stddev = 1. +space flexibility mean = None +space flexibility stddev = None +space asperity mean = None +space asperity stddev = None + +# freq mode = integrated wiener process +freq mode = cfm +freq start = 1e9 +freq pixel size = 1e8 +freq npix = 10 +freq fluctuations mean = 0.2 +freq fluctuations stddev = 0.2 +freq loglogavgslope mean = -4 +freq loglogavgslope stddev = 0.5 +freq flexibility mean = None +freq flexibility stddev = None +freq asperity mean = None +freq asperity stddev = None diff --git a/test/test_sky_models.py b/test/test_sky_models.py new file mode 100644 index 00000000..da2fcc41 --- /dev/null +++ b/test/test_sky_models.py @@ -0,0 +1,18 @@ +# SPDX-License-Identifier: GPL-3.0-or-later +# Copyright(C) 2020-2021 Max-Planck-Society +# Author: Philipp Arras + +import configparser + +import numpy as np +import pytest +import resolve as rve + +pmp = pytest.mark.parametrize +np.seterr(all="raise") + + +def test_build_multi_frequency_skymodel(): + cfg = configparser.ConfigParser() + cfg.read("test/mf_sky.cfg") + rve.multi_frequency_sky(cfg["sky"]) -- GitLab From 78b231cdcad5aa8567002c2ba686f9ff39deb162 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 25 Nov 2021 16:50:38 +0100 Subject: [PATCH 74/89] Add more general response --- resolve/response.py | 136 +------------------------- resolve/response_new.py | 212 ++++++++++++++++++++++++++++++++++++++++ resolve/util.py | 28 ++++++ test/test_general.py | 6 +- 4 files changed, 245 insertions(+), 137 deletions(-) create mode 100644 resolve/response_new.py diff --git a/resolve/response.py b/resolve/response.py index 14e0caff..1f743bdc 100644 --- a/resolve/response.py +++ b/resolve/response.py @@ -4,17 +4,13 @@ from functools import reduce from operator import add -from itertools import product - -import numpy as np -from ducc0.wgridder.experimental import dirty2vis, vis2dirty import nifty8 as ift +import numpy as np -from .constants import SPEEDOFLIGHT -from .global_config import epsilon, nthreads, wgridding, verbosity, np_dtype -from .irg_space import IRGSpace from .data.observation import Observation +from .irg_space import IRGSpace +from .response_new import SingleResponse from .util import my_assert, my_assert_isinstance, my_asserteq @@ -225,129 +221,3 @@ class ResponseDistributor(ift.LinearOperator): else: res = res + new return res - - -class FullResponse(ift.LinearOperator): - def __init__(self, observation, sky_domain): - raise NotImplementedError - - -class SingleResponse(ift.LinearOperator): - def __init__(self, domain, uvw, freq, mask=None, facets=(1, 1)): - my_assert_isinstance(facets, tuple) - for ii in range(2): - if domain.shape[ii] % facets[ii] != 0: - raise ValueError("nfacets needs to be divisor of npix.") - self._domain = ift.DomainTuple.make(domain) - tgt = ift.UnstructuredDomain(uvw.shape[0]), ift.UnstructuredDomain(freq.size) - self._target = ift.makeDomain(tgt) - self._capability = self.TIMES | self.ADJOINT_TIMES - self._args = { - "uvw": uvw, - "freq": freq, - "pixsize_x": self._domain[0].distances[0], - "pixsize_y": self._domain[0].distances[1], - "epsilon": epsilon(), - "do_wgridding": wgridding(), - "nthreads": nthreads(), - "flip_v": True, - } - if mask is not None: - self._args["mask"] = mask.astype(np.uint8) - self._vol = self._domain[0].scalar_dvol - self._target_dtype = np_dtype(True) - self._domain_dtype = np_dtype(False) - self._ofac = None - self._facets = facets - - def apply(self, x, mode): - self._check_input(x, mode) - one_facet = self._facets == (1, 1) - x = x.val - if mode == self.TIMES: - x = x.astype(self._domain_dtype, copy=False) - if one_facet: - res = self._times(x) - else: - res = self._facet_times(x) - else: - x = x.astype(self._target_dtype, copy=False) - if one_facet: - res = self._adjoint(x) - else: - res = self._facet_adjoint(x) - res = ift.makeField(self._tgt(mode), res * self._vol) - - expected_dtype = self._target_dtype if mode == self.TIMES else self._domain_dtype - my_asserteq(res.dtype, expected_dtype) - - return res - - def oversampling_factors(self): - if self._ofac is not None: - return self._ofac - maxuv = ( - np.max(np.abs(self._args["uvw"][:, 0:2]), axis=0) - * np.max(self._args["freq"]) - / SPEEDOFLIGHT - ) - hspace = self._domain[0].get_default_codomain() - hvol = np.array(hspace.shape) * np.array(hspace.distances) / 2 - self._ofac = hvol / maxuv - return self._ofac - - def _times(self, x): - if verbosity(): - print(f"\nINFO: Oversampling factors in response: {self.oversampling_factors()}\n") - return dirty2vis(dirty=x, verbosity=verbosity(), **self._args) - - def _adjoint(self, x): - my_assert(x.flags["C_CONTIGUOUS"]) - nx, ny = self._domain.shape - if verbosity(): - print(f"\nINFO: Oversampling factors in response: {self.oversampling_factors()}\n") - return vis2dirty(vis=x, npix_x=nx, npix_y=ny, verbosity=verbosity(), **self._args) - - def _facet_times(self, x): - nfacets_x, nfacets_y = self._facets - npix_x, npix_y = self._domain.shape - nx = npix_x // nfacets_x - ny = npix_y // nfacets_y - vis = None - for xx, yy in product(range(nfacets_x), range(nfacets_y)): - cx = ((0.5 + xx) / nfacets_x - 0.5) * self._args["pixsize_x"] * npix_x - cy = ((0.5 + yy) / nfacets_y - 0.5) * self._args["pixsize_y"] * npix_y - facet = x[nx * xx: nx * (xx + 1), ny * yy: ny * (yy + 1)] - foo = dirty2vis( - dirty=facet, - center_x=cx, - center_y=cy, - verbosity=verbosity(), - **self._args - ) - if vis is None: - vis = foo - else: - vis += foo - return vis - - def _facet_adjoint(self, x): - nfacets_x, nfacets_y = self._facets - npix_x, npix_y = self._domain.shape - nx = npix_x // nfacets_x - ny = npix_y // nfacets_y - res = np.zeros((npix_x, npix_y), self._domain_dtype) - for xx, yy in product(range(nfacets_x), range(nfacets_y)): - cx = ((0.5 + xx) / nfacets_x - 0.5) * self._args["pixsize_x"] * npix_x - cy = ((0.5 + yy) / nfacets_y - 0.5) * self._args["pixsize_y"] * npix_y - im = vis2dirty( - vis=x, - npix_x=nx, - npix_y=ny, - center_x=cx, - center_y=cy, - verbosity=verbosity(), - **self._args - ) - res[nx * xx: nx * (xx + 1), ny * yy: ny * (yy + 1)] = im - return res diff --git a/resolve/response_new.py b/resolve/response_new.py new file mode 100644 index 00000000..09f21e6f --- /dev/null +++ b/resolve/response_new.py @@ -0,0 +1,212 @@ +# SPDX-License-Identifier: GPL-3.0-or-later +# Copyright(C) 2019-2021 Max-Planck-Society +# Author: Philipp Arras, Jakob Knollmüller + +from itertools import product + +import nifty8 as ift +import numpy as np +from ducc0.wgridder.experimental import dirty2vis, vis2dirty + +from .constants import SPEEDOFLIGHT +from .data.observation import Observation +from .global_config import epsilon, np_dtype, nthreads, verbosity, wgridding +from .util import (assert_sky_domain, my_assert, my_assert_isinstance, + my_asserteq) + + +class InterferometryResponse(ift.LinearOperator): + def __init__(self, observation, domain): + assert isinstance(observation, Observation) + domain = ift.DomainTuple.make(domain) + pdom, tdom, fdom, sdom = domain + + assert_sky_domain(domain) + + n_time_bins = tdom.shape[0] + n_freq_bins = fdom.shape[0] + + t_binbounds = tdom.binbounds() + f_binbounds = fdom.binbounds() + + sr = [] + row_indices, freq_indices = [], [] + for tt in range(n_time_bins): + sr_tmp, t_tmp, f_tmp = [], [], [] + oo, tind = observation.restrict_by_time(t_binbounds[tt], t_binbounds[tt+1], True) + for ff in range(n_freq_bins): + ooo, find = oo.restrict_by_freq(f_binbounds[ff], f_binbounds[ff+1], True) + if any(np.array(ooo.vis.shape) == 0): + rrr = None + else: + rrr = SingleResponse(domain[3], ooo.uvw, ooo.freq) # FIXME Future: Include mask here? + sr_tmp.append(rrr) + t_tmp.append(tind) + f_tmp.append(find) + sr.append(sr_tmp) + row_indices.append(t_tmp) + freq_indices.append(f_tmp) + self._sr = sr + self._row_indices = row_indices + self._freq_indices = freq_indices + + self._domain = domain + self._target = ift.makeDomain((domain[0],) + observation.vis.domain[1:]) + self._capability = self.TIMES | self.ADJOINT_TIMES + + def apply(self, x, mode): + self._check_input(x, mode) + x = x.val + if mode == self.TIMES: + res = self._times(x) + else: + res = self._adj_times(x) + return ift.makeField(self._tgt(mode), res) + + def _times(self, x): + res = np.empty(self.target.shape, np.complex128) + res[()] = np.nan # TEMP + for pp in range(self.domain.shape[0]): + for tt in range(self.domain.shape[1]): + for ff in range(self.domain.shape[2]): + op = self._sr[tt][ff] + if op is None: + continue + inp = ift.makeField(op.domain, x[pp, tt, ff]) + res[pp, self._row_indices[tt][ff], self._freq_indices[tt][ff]] = op(inp).val + assert np.sum(np.isnan(res)) == 0 # TEMP + return res + + def _adj_times(self, x): + res = np.zeros(self.domain.shape, np.float64) + for pp in range(self.domain.shape[0]): + for tt in range(self.domain.shape[1]): + for ff in range(self.domain.shape[2]): + op = self._sr[tt][ff] + if op is None: + continue + inp = x[pp, self._row_indices[tt][ff], self._freq_indices[tt][ff]] + inp = ift.makeField(op.target, inp) + res[pp, tt, ff] = op.adjoint(inp).val + return res + + +class SingleResponse(ift.LinearOperator): + def __init__(self, domain, uvw, freq, mask=None, facets=(1, 1)): + my_assert_isinstance(facets, tuple) + for ii in range(2): + if domain.shape[ii] % facets[ii] != 0: + raise ValueError("nfacets needs to be divisor of npix.") + self._domain = ift.DomainTuple.make(domain) + tgt = ift.UnstructuredDomain(uvw.shape[0]), ift.UnstructuredDomain(freq.size) + self._target = ift.makeDomain(tgt) + self._capability = self.TIMES | self.ADJOINT_TIMES + self._args = { + "uvw": uvw, + "freq": freq, + "pixsize_x": self._domain[0].distances[0], + "pixsize_y": self._domain[0].distances[1], + "epsilon": epsilon(), + "do_wgridding": wgridding(), + "nthreads": nthreads(), + "flip_v": True, + } + if mask is not None: + self._args["mask"] = mask.astype(np.uint8) + self._vol = self._domain[0].scalar_dvol + self._target_dtype = np_dtype(True) + self._domain_dtype = np_dtype(False) + self._ofac = None + self._facets = facets + + def apply(self, x, mode): + self._check_input(x, mode) + one_facet = self._facets == (1, 1) + x = x.val + if mode == self.TIMES: + x = x.astype(self._domain_dtype, copy=False) + if one_facet: + res = self._times(x) + else: + res = self._facet_times(x) + else: + x = x.astype(self._target_dtype, copy=False) + if one_facet: + res = self._adjoint(x) + else: + res = self._facet_adjoint(x) + res = ift.makeField(self._tgt(mode), res * self._vol) + + expected_dtype = self._target_dtype if mode == self.TIMES else self._domain_dtype + my_asserteq(res.dtype, expected_dtype) + + return res + + def oversampling_factors(self): + if self._ofac is not None: + return self._ofac + maxuv = ( + np.max(np.abs(self._args["uvw"][:, 0:2]), axis=0) + * np.max(self._args["freq"]) + / SPEEDOFLIGHT + ) + hspace = self._domain[0].get_default_codomain() + hvol = np.array(hspace.shape) * np.array(hspace.distances) / 2 + self._ofac = hvol / maxuv + return self._ofac + + def _times(self, x): + if verbosity(): + print(f"\nINFO: Oversampling factors in response: {self.oversampling_factors()}\n") + return dirty2vis(dirty=x, verbosity=verbosity(), **self._args) + + def _adjoint(self, x): + my_assert(x.flags["C_CONTIGUOUS"]) + nx, ny = self._domain.shape + if verbosity(): + print(f"\nINFO: Oversampling factors in response: {self.oversampling_factors()}\n") + return vis2dirty(vis=x, npix_x=nx, npix_y=ny, verbosity=verbosity(), **self._args) + + def _facet_times(self, x): + nfacets_x, nfacets_y = self._facets + npix_x, npix_y = self._domain.shape + nx = npix_x // nfacets_x + ny = npix_y // nfacets_y + vis = None + for xx, yy in product(range(nfacets_x), range(nfacets_y)): + cx = ((0.5 + xx) / nfacets_x - 0.5) * self._args["pixsize_x"] * npix_x + cy = ((0.5 + yy) / nfacets_y - 0.5) * self._args["pixsize_y"] * npix_y + facet = x[nx * xx: nx * (xx + 1), ny * yy: ny * (yy + 1)] + foo = dirty2vis( + dirty=facet, + center_x=cx, + center_y=cy, + verbosity=verbosity(), + **self._args + ) + if vis is None: + vis = foo + else: + vis += foo + return vis + + def _facet_adjoint(self, x): + nfacets_x, nfacets_y = self._facets + npix_x, npix_y = self._domain.shape + nx = npix_x // nfacets_x + ny = npix_y // nfacets_y + res = np.zeros((npix_x, npix_y), self._domain_dtype) + for xx, yy in product(range(nfacets_x), range(nfacets_y)): + cx = ((0.5 + xx) / nfacets_x - 0.5) * self._args["pixsize_x"] * npix_x + cy = ((0.5 + yy) / nfacets_y - 0.5) * self._args["pixsize_y"] * npix_y + im = vis2dirty( + vis=x, + npix_x=nx, + npix_y=ny, + center_x=cx, + center_y=cy, + verbosity=verbosity(), + **self._args + ) + res[nx * xx: nx * (xx + 1), ny * yy: ny * (yy + 1)] = im + return res diff --git a/resolve/util.py b/resolve/util.py index 9f015fdc..d603f865 100644 --- a/resolve/util.py +++ b/resolve/util.py @@ -154,3 +154,31 @@ and `sdom` is a two-dimensional `RGSpace`. my_assert_isinstance(fdom, IRGSpace) my_assert_isinstance(sdom, ift.RGSpace) my_asserteq(len(sdom.shape), 2) + + +def assert_data_domain(dom): + """Check that input fulfils resolve's conventions of a data domain. + + A data domain is a DomainTuple: + + dom = (pdom, rdom, fdom) + + where `pdom` is a `PolarizationSpace`, `rdom` is an UnstructuredDomain representing the rows of the measurement set and + `fdom` is an `IRGSpace` containing the frequencies. + + Parameters + ---------- + dom : DomainTuple + Object that is checked to fulfil the properties. + """ + from .irg_space import IRGSpace + from .polarization_space import PolarizationSpace + + my_assert_isinstance(dom, ift.DomainTuple) + my_asserteq(len(dom), 3) + pdom, rdom, fdom = dom + for ii in range(3): + my_asserteq(len(dom[ii].shape), 1) + my_assert_isinstance(pdom, PolarizationSpace) + my_assert_isinstance(rdom, ift.UnstructuredDomain) + my_assert_isinstance(fdom, IRGSpace) diff --git a/test/test_general.py b/test/test_general.py index 3405b254..76953c32 100644 --- a/test/test_general.py +++ b/test/test_general.py @@ -274,8 +274,7 @@ def test_response_distributor(): @pmp("facets", FACETS) def test_single_response(obs, facets): mask = obs.mask.val - op = rve.response.SingleResponse(dom, obs.uvw, obs.freq, mask[0], - facets=facets) + op = rve.SingleResponse(dom, obs.uvw, obs.freq, mask[0], facets=facets) ift.extra.check_linear_operator(op, np.float64, np.complex64, only_r_linear=True, rtol=1e-6, atol=1e-6) @@ -285,8 +284,7 @@ def test_facet_consistency(): res0 = None pos = ift.from_random(dom) for facets in FACETS: - op = rve.response.SingleResponse(dom, obs.uvw, obs.freq, obs.mask.val[0], - facets=facets) + op = rve.SingleResponse(dom, obs.uvw, obs.freq, obs.mask.val[0], facets=facets) res = op(pos) if res0 is None: res0 = res -- GitLab From 39983be383ee8527db71bcc2cded01936a3e739c Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Thu, 25 Nov 2021 16:52:03 +0100 Subject: [PATCH 75/89] Work on new sky model interface --- demo/imaging_with_automatic_weighting.py | 6 ++--- resolve/__init__.py | 29 ++++++++++++------------ resolve/sky_model.py | 15 ++++++++++-- 3 files changed, 30 insertions(+), 20 deletions(-) diff --git a/demo/imaging_with_automatic_weighting.py b/demo/imaging_with_automatic_weighting.py index 5fcf278f..eaecdb8b 100644 --- a/demo/imaging_with_automatic_weighting.py +++ b/demo/imaging_with_automatic_weighting.py @@ -42,13 +42,13 @@ def main(cfg_file_name): # Sky model sky, operators = rve.single_frequency_sky(cfg["sky"]) + raw_sky = operators["sky"] enable_points = operators["points"] is not None - operators["sky"] = sky - operators["logsky"] = sky.log() + operators["logsky"] = raw_sky.log() p = ift.Plot() for _ in range(9): - p.add(sky(ift.from_random(sky.domain)), norm=LogNorm()) + p.add(raw_sky(ift.from_random(raw_sky.domain)), norm=LogNorm()) if master: p.output(name="sky_prior_samples.png") # /Sky model diff --git a/resolve/__init__.py b/resolve/__init__.py index 9613ce36..cebc6dc5 100644 --- a/resolve/__init__.py +++ b/resolve/__init__.py @@ -1,29 +1,28 @@ -from .data.antenna_positions import AntennaPositions from .calibration import calibration_distribution -from .library.calibrators import * -from .library.primary_beams import * from .constants import * -from .data.direction import * +from .data.antenna_positions import AntennaPositions from .data.averaging import * +from .data.direction import * +from .data.ms_import import * +from .data.observation import * +from .data.polarization import Polarization +from .extra import mpi_load from .fits import field2fits, fits2field from .global_config import * +from .integrated_wiener_process import (IntWProcessInitialConditions, + WienerIntegrations) +from .irg_space import IRGSpace +from .library.calibrators import * +from .library.primary_beams import * from .likelihood import * from .mosaicing import * from .mpi import onlymaster from .mpi_operators import * -from .data.ms_import import * -from .irg_space import IRGSpace -from .integrated_wiener_process import ( - IntWProcessInitialConditions, - WienerIntegrations, -) -from .data.observation import * from .points import PointInserter -from .data.polarization import Polarization from .polarization_matrix_exponential import * from .polarization_space import * -from .response import MfResponse, ResponseDistributor, StokesIResponse, SingleResponse +from .response import MfResponse, ResponseDistributor, StokesIResponse +from .response_new import InterferometryResponse, SingleResponse from .simple_operators import * -from .util import * -from .extra import mpi_load from .sky_model import * +from .util import * diff --git a/resolve/sky_model.py b/resolve/sky_model.py index 5751a718..03d127bc 100644 --- a/resolve/sky_model.py +++ b/resolve/sky_model.py @@ -40,6 +40,10 @@ def single_frequency_sky(cfg_section): add["points"] = None else: raise ValueError(f"In order to disable point source component, set `point sources mode` to `disable`. Got: {mode}") + + add["sky"] = sky + sky = DomainChangerAndReshaper(sky.target, _default_sky_domain(sdom=dom)) @ sky + assert_sky_domain(sky.target) return sky, add @@ -60,8 +64,10 @@ def multi_frequency_sky(cfg): add["freq normalized power spectrum"] = fampl**2 add["space normalized power spectrum"] = sampl**2 - tgt = PolarizationSpace("I"), IRGSpace([np.nan]), fdom, sdom - sky = DomainChangerAndReshaper(sky.target, tgt) @ sky + add["sky"] = sky + reshape = DomainChangerAndReshaper(sky.target, _default_sky_domain(fdom=fdom, sdom=sdom)) + sky = reshape @ sky + elif cfg["freq mode"] == "integrated wiener process": raise NotImplementedError else: @@ -111,3 +117,8 @@ def _parse_or_none(cfg, key): return None if key0 in cfg: return (cfg.getfloat(key0), cfg.getfloat(key1)) + + +def _default_sky_domain(pdom=PolarizationSpace(), tdom=IRGSpace([np.nan]), fdom=IRGSpace([np.nan]), + sdom=ift.RGSpace([1, 1])): + return pdom, tdom, fdom, sdom -- GitLab From a426868546774bc98d18c2eca31c31c6e6d7ebff Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 26 Nov 2021 15:10:23 +0100 Subject: [PATCH 76/89] Fixup --- resolve/sky_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/resolve/sky_model.py b/resolve/sky_model.py index 03d127bc..1d760e90 100644 --- a/resolve/sky_model.py +++ b/resolve/sky_model.py @@ -119,6 +119,6 @@ def _parse_or_none(cfg, key): return (cfg.getfloat(key0), cfg.getfloat(key1)) -def _default_sky_domain(pdom=PolarizationSpace(), tdom=IRGSpace([np.nan]), fdom=IRGSpace([np.nan]), +def _default_sky_domain(pdom=PolarizationSpace("I"), tdom=IRGSpace([np.nan]), fdom=IRGSpace([np.nan]), sdom=ift.RGSpace([1, 1])): return pdom, tdom, fdom, sdom -- GitLab From ac97d79411c0f5a57c6e811b6453795b622efd9f Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 26 Nov 2021 15:13:59 +0100 Subject: [PATCH 77/89] Reshaper, DomainChanger -> DomainChangerAndReshaper --- demo/calibration_and_imaging.py | 2 +- resolve/util.py | 21 --------------------- test/test_general.py | 2 +- 3 files changed, 2 insertions(+), 23 deletions(-) diff --git a/demo/calibration_and_imaging.py b/demo/calibration_and_imaging.py index 319500ce..e8234cb6 100644 --- a/demo/calibration_and_imaging.py +++ b/demo/calibration_and_imaging.py @@ -33,7 +33,7 @@ def main(): good_size(int(2 * (tmax - tmin) / solution_interval)), solution_interval ) print(f"Npix in time domain {time_domain.shape[0]}") - reshaper = rve.Reshaper( + reshaper = rve.DomainChangerAndReshaper( [ift.UnstructuredDomain(total_N), time_domain], [ ift.UnstructuredDomain(npol), diff --git a/resolve/util.py b/resolve/util.py index d603f865..828c0994 100644 --- a/resolve/util.py +++ b/resolve/util.py @@ -69,17 +69,6 @@ def _types_equal(a, b): return type(a) == type(b) -class Reshaper(ift.LinearOperator): - def __init__(self, domain, target): - self._domain = ift.DomainTuple.make(domain) - self._target = ift.DomainTuple.make(target) - self._capability = self.TIMES | self.ADJOINT_TIMES - - def apply(self, x, mode): - self._check_input(x, mode) - return ift.makeField(self._tgt(mode), x.val.reshape(self._tgt(mode).shape)) - - def divide_where_possible(a, b): # Otherwise one if isinstance(a, ift.Field): @@ -115,16 +104,6 @@ def rows_to_baselines(antenna_positions, data_field): return res -class DomainChanger(ift.LinearOperator): - def __init__(self, domain, target): - self._domain = ift.makeDomain(domain) - self._target = ift.makeDomain(target) - self._capability = self.TIMES | self.ADJOINT_TIMES - - def apply(self, x, mode): - return ift.makeField(self._tgt(mode), x.val) - - def assert_sky_domain(dom): """Check that input fulfils resolve's conventions of a sky domain. diff --git a/test/test_general.py b/test/test_general.py index 76953c32..78e31151 100644 --- a/test/test_general.py +++ b/test/test_general.py @@ -161,7 +161,7 @@ def test_calibration_likelihood(time_mode): ift.UnstructuredDomain(nfreq), ] if time_mode: - reshaper = rve.Reshaper([ift.UnstructuredDomain(total_N), time_domain], dom) + reshaper = rve.DomainChangerAndReshaper([ift.UnstructuredDomain(total_N), time_domain], dom) dct = {"offset_mean": 0, "offset_std": (1, 0.5)} dct1 = { "fluctuations": (2.0, 1.0), -- GitLab From 640609c6f642fe354118b6167c398f35cd7a472a Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 26 Nov 2021 16:19:07 +0100 Subject: [PATCH 78/89] Add wiener process model --- resolve/sky_model.py | 70 +++++++++++++++++++++++++++++++++- test/mf_sky.cfg | 4 +- test/mf_sky_wiener_process.cfg | 41 ++++++++++++++++++++ test/test_sky_models.py | 13 ++++++- 4 files changed, 123 insertions(+), 5 deletions(-) create mode 100644 test/mf_sky_wiener_process.cfg diff --git a/resolve/sky_model.py b/resolve/sky_model.py index 1d760e90..c0b705a6 100644 --- a/resolve/sky_model.py +++ b/resolve/sky_model.py @@ -6,6 +6,8 @@ import nifty8 as ift import numpy as np from .constants import str2rad +from .integrated_wiener_process import (IntWProcessInitialConditions, + WienerIntegrations) from .irg_space import IRGSpace from .points import PointInserter from .polarization_space import PolarizationSpace @@ -52,6 +54,7 @@ def multi_frequency_sky(cfg): raise NotImplementedError("Point sources are not supported yet.") sdom = _spatial_dom(cfg) + add = {} if cfg["freq mode"] == "cfm": fnpix, df = cfg.getfloat("freq npix"), cfg.getfloat("freq pixel size") @@ -59,7 +62,7 @@ def multi_frequency_sky(cfg): fdom_rg = ift.RGSpace(fnpix, df) logdiffuse, cfm = cfm_from_cfg(cfg, {"freq": fdom_rg, "space": sdom}, "sky diffuse") sky = logdiffuse.exp() - add = {"logdiffuse": logdiffuse} + add["logdiffuse"] = logdiffuse fampl, sampl = list(cfm.get_normalized_amplitudes()) add["freq normalized power spectrum"] = fampl**2 add["space normalized power spectrum"] = sampl**2 @@ -69,7 +72,70 @@ def multi_frequency_sky(cfg): sky = reshape @ sky elif cfg["freq mode"] == "integrated wiener process": - raise NotImplementedError + prefix = "sky" + + freq = np.sort(np.array(list(map(float, cfg["frequencies"].split(","))))) + fdom = IRGSpace(freq) + log_fdom = IRGSpace(np.log(freq / freq.mean())) + + # Base sky brightness distribution for mean frequency \nu_0 + i_0, i_0_cfm = cfm_from_cfg(cfg, {"space i0": sdom}, "space i0") + + # simple spectral index that is frequency independent + alpha, alpha_cfm = cfm_from_cfg(cfg, {"space alpha": sdom}, "space alpha") + + flexibility = _parse_or_none(cfg, "freq flexibility") + if flexibility is None: + raise RuntimeError("freq flexibility cannot be None") + asperity = _parse_or_none(cfg, "freq asperity") + + # IDEA Try to use individual power spectra + n_freq_xi_fields = 2 * (log_fdom.size - 1) + cfm = ift.CorrelatedFieldMaker(f"{prefix}freq_xi", total_N=n_freq_xi_fields) + cfm.set_amplitude_total_offset(0.0, None) + # FIXME Support fixed fluctuations + cfm.add_fluctuations(sdom, (1, 1e-6), (1.2, 0.4), (0.2, 0.2), (-2, 0.5), dofdex=n_freq_xi_fields * [0]) + freq_xi = cfm.finalize(0) + + # Integrate over excitation fields + intop = WienerIntegrations(log_fdom, sdom) + freq_xi = DomainChangerAndReshaper(freq_xi.target, intop.domain) @ freq_xi + broadcast = ift.ContractionOperator(intop.domain[0], None).adjoint + broadcast_full = ift.ContractionOperator(intop.domain, 1).adjoint + vol = log_fdom.distances + + flex = ift.LognormalTransform(*flexibility, prefix + "flexibility", 0) + + dom = intop.domain[0] + vflex = np.empty(dom.shape) + vflex[0] = vflex[1] = np.sqrt(vol) + sig_flex = ift.makeOp(ift.makeField(dom, vflex)) @ broadcast @ flex + sig_flex = broadcast_full @ sig_flex + shift = np.ones(dom.shape) + shift[0] = vol * vol / 12.0 + if asperity is None: + shift = ift.DiagonalOperator(ift.makeField(dom, shift).sqrt(), intop.domain, 0) + increments = shift @ (freq_xi * sig_flex) + else: + asp = ift.LognormalTransform(*asperity, prefix + "asperity", 0) + vasp = np.empty(dom.shape) + vasp[0] = 1 + vasp[1] = 0 + vasp = ift.DiagonalOperator(ift.makeField(dom, vasp), domain=broadcast.target, spaces=0) + sig_asp = broadcast_full @ vasp @ broadcast @ asp + shift = ift.makeField(intop.domain, np.broadcast_to(shift[..., None, None], intop.domain.shape)) + increments = freq_xi * sig_flex * (ift.Adder(shift) @ sig_asp).ptw("sqrt") + logsky = IntWProcessInitialConditions(i_0, alpha, intop @ increments) + + sky = logsky.exp() + add["sky"] = sky + add["logsky"] = logsky + add["i0"] = i_0 + add["alpha"] = alpha + + reshape = DomainChangerAndReshaper(sky.target, _default_sky_domain(fdom=fdom, sdom=sdom)) + sky = reshape @ sky + else: raise RuntimeError diff --git a/test/mf_sky.cfg b/test/mf_sky.cfg index d486869e..32b0558d 100644 --- a/test/mf_sky.cfg +++ b/test/mf_sky.cfg @@ -1,8 +1,8 @@ [sky] space npix x = 500 space npix y = 500 -space fov x = 1000muas -space fov y = 1000muas +space fov x = 1deg +space fov y = 1deg # point sources mode = single # point sources locations = 0deg,0deg;0.35as,-0.22as diff --git a/test/mf_sky_wiener_process.cfg b/test/mf_sky_wiener_process.cfg new file mode 100644 index 00000000..521546af --- /dev/null +++ b/test/mf_sky_wiener_process.cfg @@ -0,0 +1,41 @@ +[sky] +space npix x = 500 +space npix y = 500 +space fov x = 1deg +space fov y = 1deg + +# point sources mode = single +# point sources locations = 0deg,0deg;0.35as,-0.22as +# point sources alpha = 0.5 +# point sources q = 0.2 + +space i0 zero mode offset = 0 +space i0 zero mode mean = 1 +space i0 zero mode stddev = 0.3 +space i0 fluctuations mean = 1.5 +space i0 fluctuations stddev = 1. +space i0 loglogavgslope mean = -3 +space i0 loglogavgslope stddev = 1. +space i0 flexibility mean = None +space i0 flexibility stddev = None +space i0 asperity mean = None +space i0 asperity stddev = None + +space alpha zero mode offset = 0 +space alpha zero mode mean = 1 +space alpha zero mode stddev = 0.3 +space alpha fluctuations mean = 1.5 +space alpha fluctuations stddev = 1. +space alpha loglogavgslope mean = -3 +space alpha loglogavgslope stddev = 1. +space alpha flexibility mean = None +space alpha flexibility stddev = None +space alpha asperity mean = None +space alpha asperity stddev = None + +freq mode = integrated wiener process +frequencies = 1000e8,1003e8,1100e8,1324e8 +freq flexibility mean = 1 +freq flexibility stddev = 1 +freq asperity mean = None +freq asperity stddev = None diff --git a/test/test_sky_models.py b/test/test_sky_models.py index da2fcc41..10d8f948 100644 --- a/test/test_sky_models.py +++ b/test/test_sky_models.py @@ -5,6 +5,7 @@ import configparser import numpy as np +import nifty8 as ift import pytest import resolve as rve @@ -15,4 +16,14 @@ np.seterr(all="raise") def test_build_multi_frequency_skymodel(): cfg = configparser.ConfigParser() cfg.read("test/mf_sky.cfg") - rve.multi_frequency_sky(cfg["sky"]) + op, _ = rve.multi_frequency_sky(cfg["sky"]) + op(ift.from_random(op.domain)) + + cfg.read("test/mf_sky_wiener_process.cfg") + op, _ = rve.multi_frequency_sky(cfg["sky"]) + op(ift.from_random(op.domain)) + + cfg["sky"]["freq asperity mean"] = "0.1" + cfg["sky"]["freq asperity stddev"] = "0.1" + op, _ = rve.multi_frequency_sky(cfg["sky"]) + op(ift.from_random(op.domain)) -- GitLab From ab3eb307da76c04aeb043312670537aaaf61814f Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 26 Nov 2021 17:02:17 +0100 Subject: [PATCH 79/89] Response: introduce shortcuts --- resolve/response_new.py | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/resolve/response_new.py b/resolve/response_new.py index 09f21e6f..06a6eb4e 100644 --- a/resolve/response_new.py +++ b/resolve/response_new.py @@ -19,10 +19,10 @@ class InterferometryResponse(ift.LinearOperator): def __init__(self, observation, domain): assert isinstance(observation, Observation) domain = ift.DomainTuple.make(domain) - pdom, tdom, fdom, sdom = domain - assert_sky_domain(domain) + pdom, tdom, fdom, sdom = domain + n_time_bins = tdom.shape[0] n_freq_bins = fdom.shape[0] @@ -33,7 +33,11 @@ class InterferometryResponse(ift.LinearOperator): row_indices, freq_indices = [], [] for tt in range(n_time_bins): sr_tmp, t_tmp, f_tmp = [], [], [] - oo, tind = observation.restrict_by_time(t_binbounds[tt], t_binbounds[tt+1], True) + if tuple(t_binbounds[tt:tt+2]) == (-np.inf, np.inf): + oo = observation + tind = slice(None) + else: + oo, tind = observation.restrict_by_time(t_binbounds[tt], t_binbounds[tt+1], True) for ff in range(n_freq_bins): ooo, find = oo.restrict_by_freq(f_binbounds[ff], f_binbounds[ff+1], True) if any(np.array(ooo.vis.shape) == 0): @@ -54,6 +58,14 @@ class InterferometryResponse(ift.LinearOperator): self._target = ift.makeDomain((domain[0],) + observation.vis.domain[1:]) self._capability = self.TIMES | self.ADJOINT_TIMES + foo = np.zeros(self.target.shape, np.int8) + for pp in range(self.domain.shape[0]): + for tt in range(self.domain.shape[1]): + for ff in range(self.domain.shape[2]): + foo[pp, self._row_indices[tt][ff], self._freq_indices[tt][ff]] = 1. + if np.any(foo == 0): + raise RuntimeError("This should not happen. This is a bug in `InterferometryResponse`.") + def apply(self, x, mode): self._check_input(x, mode) x = x.val @@ -65,7 +77,6 @@ class InterferometryResponse(ift.LinearOperator): def _times(self, x): res = np.empty(self.target.shape, np.complex128) - res[()] = np.nan # TEMP for pp in range(self.domain.shape[0]): for tt in range(self.domain.shape[1]): for ff in range(self.domain.shape[2]): @@ -74,7 +85,6 @@ class InterferometryResponse(ift.LinearOperator): continue inp = ift.makeField(op.domain, x[pp, tt, ff]) res[pp, self._row_indices[tt][ff], self._freq_indices[tt][ff]] = op(inp).val - assert np.sum(np.isnan(res)) == 0 # TEMP return res def _adj_times(self, x): -- GitLab From b2d27c6897e1fb5d03c4e50a9b2fee117ca9577b Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 26 Nov 2021 17:02:38 +0100 Subject: [PATCH 80/89] Likelihood: Adopt to new response interface --- resolve/likelihood.py | 52 +++++++++---------------------------------- 1 file changed, 11 insertions(+), 41 deletions(-) diff --git a/resolve/likelihood.py b/resolve/likelihood.py index 642c7159..11892aa1 100644 --- a/resolve/likelihood.py +++ b/resolve/likelihood.py @@ -106,11 +106,7 @@ def ImagingLikelihood( inverse_covariance_operator=None, calibration_operator=None, ): - """Versatile likelihood class that automatically chooses the correct - response class. - - Supports polarization imaging and Stokes I imaging. Supports single-frequency - and multi-frequency imaging. + """Versatile likelihood class. If a calibration operator is passed, it returns an operator that computes: @@ -133,11 +129,12 @@ def ImagingLikelihood( observation.weight is used for computing the likelihood. sky_operator : Operator - Operator that generates sky. Can have as a target either a two-dimensional - DomainTuple (single-frequency imaging), a three-dimensional DomainTuple - (multi-frequency imaging) or a MultiDomain (polarization imaging). - For multi-frequency imaging it is required that the first entry of the - three-dimensional DomainTuple represents the domain of the frequencies. + Operator that generates sky. Needs to have as target: + + dom = (pdom, tdom, fdom, sdom) + + where `pdom` is a `PolarizationSpace`, `tdom` and `fdom` are an + `IRGSpace`, and `sdom` is a two-dimensional `RGSpace`. inverse_covariance_operator : Operator Optional. Target needs to be the same space as observation.vis. If it is @@ -145,38 +142,12 @@ def ImagingLikelihood( calibration_operator : Operator Optional. Target needs to be the same as observation.vis. + """ my_assert_isinstance(sky_operator, ift.Operator) - sdom = sky_operator.target - - mosaicing = isinstance(observation, dict) - - if isinstance(sdom, ift.MultiDomain) and not mosaicing: - if len(sdom["I"].shape) == 3: - raise NotImplementedError( - "Polarization and multi-frequency at the same time not supported yet." - ) - else: - R = FullPolResponse(observation, sky_operator.target) - else: - if not mosaicing and len(sdom.shape) == 3: - R = MfResponse(observation, sdom[0], sdom[1]) - else: - R = StokesIResponse(observation, sdom) - model_data = R @ sky_operator + model_data = InterferometryResponse(observation, sky_operator.target) @ sky_operator if inverse_covariance_operator is None: - if mosaicing: - vis, weight, mask = {}, {}, {} - for kk, oo in observation.items(): - vis[kk], weight[kk] = oo.vis, oo.weight - vis = ift.MultiField.from_dict(vis) - weight = ift.MultiField.from_dict(weight) - # FIXME Use mask here as well - mask = get_mask_multi_field(weight) - vis = mask(vis) - weight = mask(weight) - else: - mask, vis, weight = _get_mask(observation) + mask, vis, weight = _get_mask(observation) return _build_gauss_lh_nres(mask @ model_data, vis, weight) return _varcov(observation, model_data, inverse_covariance_operator) @@ -187,8 +158,7 @@ def CalibrationLikelihood( model_visibilities, inverse_covariance_operator=None, ): - """Versatile calibration likelihood class that automatically chooses - the correct response class. + """Versatile calibration likelihood class It returns an operator that computes: -- GitLab From ef54388c552b04b36b7c78494324caced790490a Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 26 Nov 2021 17:03:00 +0100 Subject: [PATCH 81/89] Cosmetics --- resolve/likelihood.py | 5 ++--- resolve/response_new.py | 20 ++++---------------- 2 files changed, 6 insertions(+), 19 deletions(-) diff --git a/resolve/likelihood.py b/resolve/likelihood.py index 11892aa1..1d5028ea 100644 --- a/resolve/likelihood.py +++ b/resolve/likelihood.py @@ -5,12 +5,11 @@ from functools import reduce from operator import add -import numpy as np - import nifty8 as ift +import numpy as np from .data.observation import Observation -from .response import FullPolResponse, MfResponse, StokesIResponse +from .response_new import InterferometryResponse from .util import my_assert, my_assert_isinstance, my_asserteq diff --git a/resolve/response_new.py b/resolve/response_new.py index 06a6eb4e..5546dbbd 100644 --- a/resolve/response_new.py +++ b/resolve/response_new.py @@ -187,13 +187,8 @@ class SingleResponse(ift.LinearOperator): cx = ((0.5 + xx) / nfacets_x - 0.5) * self._args["pixsize_x"] * npix_x cy = ((0.5 + yy) / nfacets_y - 0.5) * self._args["pixsize_y"] * npix_y facet = x[nx * xx: nx * (xx + 1), ny * yy: ny * (yy + 1)] - foo = dirty2vis( - dirty=facet, - center_x=cx, - center_y=cy, - verbosity=verbosity(), - **self._args - ) + foo = dirty2vis(dirty=facet, center_x=cx, center_y=cy, verbosity=verbosity(), + **self._args) if vis is None: vis = foo else: @@ -209,14 +204,7 @@ class SingleResponse(ift.LinearOperator): for xx, yy in product(range(nfacets_x), range(nfacets_y)): cx = ((0.5 + xx) / nfacets_x - 0.5) * self._args["pixsize_x"] * npix_x cy = ((0.5 + yy) / nfacets_y - 0.5) * self._args["pixsize_y"] * npix_y - im = vis2dirty( - vis=x, - npix_x=nx, - npix_y=ny, - center_x=cx, - center_y=cy, - verbosity=verbosity(), - **self._args - ) + im = vis2dirty(vis=x, npix_x=nx, npix_y=ny, center_x=cx, center_y=cy, + verbosity=verbosity(), **self._args) res[nx * xx: nx * (xx + 1), ny * yy: ny * (yy + 1)] = im return res -- GitLab From 24d250c058e06d274c15727bc60d1c04212fe1d2 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 26 Nov 2021 17:03:44 +0100 Subject: [PATCH 82/89] Skymodel: compatibility with new likelihood interface --- resolve/sky_model.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/resolve/sky_model.py b/resolve/sky_model.py index c0b705a6..6fbc0358 100644 --- a/resolve/sky_model.py +++ b/resolve/sky_model.py @@ -44,7 +44,10 @@ def single_frequency_sky(cfg_section): raise ValueError(f"In order to disable point source component, set `point sources mode` to `disable`. Got: {mode}") add["sky"] = sky - sky = DomainChangerAndReshaper(sky.target, _default_sky_domain(sdom=dom)) @ sky + conv = DomainChangerAndReshaper(sky.target, _default_sky_domain(sdom=dom)) + sky = conv @ sky + if add["points"] is not None: + add["points"] = conv @ add["points"] assert_sky_domain(sky.target) return sky, add -- GitLab From 2aaa3519d1208fc403ee9d8020f264e2f4543d63 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 26 Nov 2021 17:04:35 +0100 Subject: [PATCH 83/89] Use double precision by default --- resolve/global_config.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/resolve/global_config.py b/resolve/global_config.py index f9d06374..ba478131 100644 --- a/resolve/global_config.py +++ b/resolve/global_config.py @@ -9,7 +9,7 @@ _wgridding = None _epsilon = None _nthreads = 1 _verbosity = False -_double_precision = False +_double_precision = True def wgridding(): -- GitLab From 2465693be954b91030719c31b7b19d195b6b57f1 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Fri, 26 Nov 2021 17:05:58 +0100 Subject: [PATCH 84/89] rve.KeyPrefixer -> ift.PrependKey --- resolve/likelihood.py | 5 ++--- resolve/simple_operators.py | 25 ++----------------------- test/test_general.py | 8 -------- 3 files changed, 4 insertions(+), 34 deletions(-) diff --git a/resolve/likelihood.py b/resolve/likelihood.py index 1d5028ea..51fb1d4c 100644 --- a/resolve/likelihood.py +++ b/resolve/likelihood.py @@ -59,7 +59,6 @@ def _build_gauss_lh_nres(op, mean, invcov): def _varcov(observation, Rs, inverse_covariance_operator): - from .simple_operators import KeyPrefixer mosaicing = isinstance(observation, dict) s0, s1 = "residual", "inverse covariance" if mosaicing: @@ -78,8 +77,8 @@ def _varcov(observation, Rs, inverse_covariance_operator): lhs.append(e @ (a+b)) masks = reduce(add, masks) - a = KeyPrefixer(masks.target, "modeld") @ masks @ Rs - b = KeyPrefixer(masks.target, "icov") @ masks @ inverse_covariance_operator + a = ift.PrependKey(masks.target, "modeld") @ masks @ Rs + b = ift.PrependKey(masks.target, "icov") @ masks @ inverse_covariance_operator lh = reduce(add, lhs) @ (a + b) vis = ift.MultiField.from_dict(vis) diff --git a/resolve/simple_operators.py b/resolve/simple_operators.py index 8be64027..4f325f3e 100644 --- a/resolve/simple_operators.py +++ b/resolve/simple_operators.py @@ -69,27 +69,6 @@ class AddEmptyDimensionAtEnd(ift.LinearOperator): return ift.makeField(self._tgt(mode), x) -class KeyPrefixer(ift.LinearOperator): - def __init__(self, domain, prefix): - self._domain = ift.MultiDomain.make(domain) - self._target = ift.MultiDomain.make( - {prefix + kk: vv for kk, vv in self._domain.items()} - ) - self._capability = self.TIMES | self.ADJOINT_TIMES - self._prefix = prefix - - def apply(self, x, mode): - self._check_input(x, mode) - if mode == self.TIMES: - res = {self._prefix + kk: vv for kk, vv in x.items()} - else: - res = {kk[len(self._prefix) :]: vv for kk, vv in x.items()} - return ift.MultiField.from_dict(res) - - def __repr__(self): - return f"{self.domain.keys()} -> {self.target.keys()}" - - def MultiDomainVariableCovarianceGaussianEnergy(data, signal_response, invcov): from .likelihood import get_mask_multi_field @@ -110,8 +89,8 @@ def MultiDomainVariableCovarianceGaussianEnergy(data, signal_response, invcov): data.domain[kk], "resi" + kk, "icov" + kk, data[kk].dtype ) ) - resi = KeyPrefixer(data.domain, "resi") @ ift.Adder(data, True) @ signal_response - invcov = KeyPrefixer(data.domain, "icov") @ invcov + resi = ift.PrependKey(data.domain, "resi") @ ift.Adder(data, True) @ signal_response + invcov = ift.PrependKey(data.domain, "icov") @ invcov return reduce(add, res) @ (resi + invcov) diff --git a/test/test_general.py b/test/test_general.py index 78e31151..a5cdb6a1 100644 --- a/test/test_general.py +++ b/test/test_general.py @@ -324,14 +324,6 @@ def test_intop(): ift.extra.check_linear_operator(op) -def test_prefixer(): - op = rve.KeyPrefixer( - ift.MultiDomain.make({"a": ift.UnstructuredDomain(10), "b": ift.RGSpace(190)}), - "invcov_inp", - ).adjoint - ift.extra.check_linear_operator(op) - - def test_global_config(): rve.set_double_precision(True) assert rve.np_dtype() == rve.np_dtype(False) == np.float64 -- GitLab From 18b1de8367d2ff0a21e74e5628e655f96c2eb15d Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Tue, 30 Nov 2021 12:25:49 +0100 Subject: [PATCH 85/89] Add calibration to likelihood --- resolve/likelihood.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/resolve/likelihood.py b/resolve/likelihood.py index 51fb1d4c..72a0fc08 100644 --- a/resolve/likelihood.py +++ b/resolve/likelihood.py @@ -144,6 +144,8 @@ def ImagingLikelihood( """ my_assert_isinstance(sky_operator, ift.Operator) model_data = InterferometryResponse(observation, sky_operator.target) @ sky_operator + if calibration_operator is not None: + model_data = calibration_operator * model_data if inverse_covariance_operator is None: mask, vis, weight = _get_mask(observation) return _build_gauss_lh_nres(mask @ model_data, vis, weight) -- GitLab From 5b3436a86bed08fd6a521cd12f40a537e4db75d5 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Tue, 30 Nov 2021 15:18:25 +0100 Subject: [PATCH 86/89] Fixups regarding polarization in likelihood --- resolve/likelihood.py | 3 +++ resolve/polarization_space.py | 23 ++++++++++++++++++++++- resolve/sky_model.py | 10 +++++++--- test/test_general.py | 5 ++++- 4 files changed, 36 insertions(+), 5 deletions(-) diff --git a/resolve/likelihood.py b/resolve/likelihood.py index 72a0fc08..aab038f1 100644 --- a/resolve/likelihood.py +++ b/resolve/likelihood.py @@ -142,8 +142,11 @@ def ImagingLikelihood( Optional. Target needs to be the same as observation.vis. """ + from .polarization_space import polarization_converter + my_assert_isinstance(sky_operator, ift.Operator) model_data = InterferometryResponse(observation, sky_operator.target) @ sky_operator + model_data = polarization_converter(model_data.target, observation.vis.domain) @ model_data if calibration_operator is not None: model_data = calibration_operator * model_data if inverse_covariance_operator is None: diff --git a/resolve/polarization_space.py b/resolve/polarization_space.py index dad3e786..6c21c883 100644 --- a/resolve/polarization_space.py +++ b/resolve/polarization_space.py @@ -20,7 +20,9 @@ class PolarizationSpace(ift.UnstructuredDomain): def __init__(self, polarization_labels): if isinstance(polarization_labels, str): - polarization_labels = (polarization_labels,) + polarization_labels = [polarization_labels] + lbl = list(polarization_labels) + lbl.sort() self._lbl = tuple(polarization_labels) for ll in self._lbl: my_assert(ll in ["I", "Q", "U", "V", "LL", "LR", "RL", "RR", "XX", "XY", "YX", "YY"]) @@ -32,3 +34,22 @@ class PolarizationSpace(ift.UnstructuredDomain): @property def labels(self): return self._lbl + + +def polarization_converter(domain, target): + from .util import my_assert_isinstance + from .simple_operators import DomainChangerAndReshaper + + domain = ift.DomainTuple.make(domain) + target = ift.DomainTuple.make(target) + my_assert_isinstance(domain[0], PolarizationSpace) + my_assert_isinstance(target[0], PolarizationSpace) + if domain is target: + return ift.ScalingOperator(domain, 1.) + + if domain[0].labels == ("I",): + if target[0].labels in [("LL", "RR"), ("XX", "YY")]: + # Convention: Stokes I 1Jy source leads to 1Jy in LL and 1Jy in RR + op = ift.ContractionOperator(target, 0).adjoint + return op @ DomainChangerAndReshaper(domain, op.domain) + raise NotImplementedError(f"Polarization converter\ndomain:\n{domain[0]}\ntarget\n{target[0]}\n") diff --git a/resolve/sky_model.py b/resolve/sky_model.py index 6fbc0358..4ab5208e 100644 --- a/resolve/sky_model.py +++ b/resolve/sky_model.py @@ -188,6 +188,10 @@ def _parse_or_none(cfg, key): return (cfg.getfloat(key0), cfg.getfloat(key1)) -def _default_sky_domain(pdom=PolarizationSpace("I"), tdom=IRGSpace([np.nan]), fdom=IRGSpace([np.nan]), - sdom=ift.RGSpace([1, 1])): - return pdom, tdom, fdom, sdom +def default_sky_domain(pdom=PolarizationSpace("I"), tdom=IRGSpace([np.nan]), fdom=IRGSpace([np.nan]), + sdom=ift.RGSpace([1, 1])): + from .util import my_assert_isinstance + + for dd in [pdom, tdom, fdom, sdom]: + my_assert_isinstance(dd, ift.Domain) + return ift.DomainTuple.make((pdom, tdom, fdom, sdom)) diff --git a/test/test_general.py b/test/test_general.py index a5cdb6a1..e795106c 100644 --- a/test/test_general.py +++ b/test/test_general.py @@ -19,7 +19,7 @@ rve.set_nthreads(nthreads) rve.set_epsilon(1e-4) rve.set_wgridding(False) OBS = [] -for polmode in ["all", "stokesi", ["LL"], "stokesiavg"]: +for polmode in ["all", "stokesi", "stokesiavg"]: OBS.append( rve.ms2observations( f"{direc}CYG-ALL-2052-2MHZ.ms", "DATA", True, 0, polarizations=polmode @@ -36,6 +36,7 @@ points = ift.InverseGammaOperator( inserter.domain, alpha=0.5, q=0.2 / dom.scalar_dvol ).ducktape("points") sky = rve.vla_beam(dom, np.mean(OBS[0].freq)) @ (sky0 + inserter @ points) +sky = rve.DomainChangerAndReshaper(sky.target, rve.default_sky_domain(sdom=sky.target[0])) @ sky freqdomain = rve.IRGSpace(np.linspace(1000, 1050, num=10)) FACETS = [(1, 1), (2, 2), (2, 1), (1, 4)] @@ -93,6 +94,8 @@ def try_lh(obs, lh_class, *args): @pmp("obs", OBS) def test_imaging_likelihood(obs): + obs = obs.restrict_to_stokesi() + op = rve.ImagingLikelihood(obs, sky) try_lh(obs, rve.ImagingLikelihood, obs, sky) -- GitLab From 9a3d95e82037714ab70b684676d3dfaa85b5a3c5 Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Tue, 30 Nov 2021 17:16:16 +0100 Subject: [PATCH 87/89] Observation: restrict_to_stokesi works on stokesi only observation --- resolve/data/observation.py | 2 ++ test/test_general.py | 1 - 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/resolve/data/observation.py b/resolve/data/observation.py index 2da44d84..3e2bcb31 100644 --- a/resolve/data/observation.py +++ b/resolve/data/observation.py @@ -490,6 +490,8 @@ class Observation(BaseObservation): return res def restrict_to_stokesi(self): + if self.vis.domain[0].labels == ("I",): + return self # FIXME Do I need to change something in self._auxiliary_tables? ind = self._polarization.stokes_i_indices() vis = self._vis[ind] diff --git a/test/test_general.py b/test/test_general.py index e795106c..54128258 100644 --- a/test/test_general.py +++ b/test/test_general.py @@ -25,7 +25,6 @@ for polmode in ["all", "stokesi", "stokesiavg"]: f"{direc}CYG-ALL-2052-2MHZ.ms", "DATA", True, 0, polarizations=polmode )[0] ) -assert OBS[1].max_snr() >= OBS[2].max_snr() # Average data so SNR increases npix, fov = 256, 1 * rve.DEG2RAD dom = ift.RGSpace((npix, npix), (fov / npix, fov / npix)) sky0 = ift.SimpleCorrelatedField( -- GitLab From 76b01d6ff3da3ce241f4d2b8e2dc7102d638deae Mon Sep 17 00:00:00 2001 From: Philipp Arras <parras@mpa-garching.mpg.de> Date: Tue, 30 Nov 2021 17:38:19 +0100 Subject: [PATCH 88/89] Fixup --- resolve/sky_model.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/resolve/sky_model.py b/resolve/sky_model.py index 4ab5208e..94013d2b 100644 --- a/resolve/sky_model.py +++ b/resolve/sky_model.py @@ -44,7 +44,7 @@ def single_frequency_sky(cfg_section): raise ValueError(f"In order to disable point source component, set `point sources mode` to `disable`. Got: {mode}") add["sky"] = sky - conv = DomainChangerAndReshaper(sky.target, _default_sky_domain(sdom=dom)) + conv = DomainChangerAndReshaper(sky.target, default_sky_domain(sdom=dom)) sky = conv @ sky if add["points"] is not None: add["points"] = conv @ add["points"] @@ -71,7 +71,7 @@ def multi_frequency_sky(cfg): add["space normalized power spectrum"] = sampl**2 add["sky"] = sky - reshape = DomainChangerAndReshaper(sky.target, _default_sky_domain(fdom=fdom, sdom=sdom)) + reshape = DomainChangerAndReshaper(sky.target, default_sky_domain(fdom=fdom, sdom=sdom)) sky = reshape @ sky elif cfg["freq mode"] == "integrated wiener process": @@ -136,7 +136,7 @@ def multi_frequency_sky(cfg): add["i0"] = i_0 add["alpha"] = alpha - reshape = DomainChangerAndReshaper(sky.target, _default_sky_domain(fdom=fdom, sdom=sdom)) + reshape = DomainChangerAndReshaper(sky.target, default_sky_domain(fdom=fdom, sdom=sdom)) sky = reshape @ sky else: -- GitLab From 53df267f0a9605f0c88c1f549fe907dbbaa55d40 Mon Sep 17 00:00:00 2001 From: Philipp Arras <c@philipp-arras.de> Date: Wed, 1 Dec 2021 17:52:22 +0100 Subject: [PATCH 89/89] Instructions for demos --- demos/README.md | 16 ++++++++++++++++ {demo => demos}/cygnusa.cfg | 0 .../imaging_with_automatic_weighting.py | 0 {demo => demos/old}/basic_polarization.py | 0 {demo => demos/old}/calibration_and_imaging.py | 0 5 files changed, 16 insertions(+) create mode 100644 demos/README.md rename {demo => demos}/cygnusa.cfg (100%) rename {demo => demos}/imaging_with_automatic_weighting.py (100%) rename {demo => demos/old}/basic_polarization.py (100%) rename {demo => demos/old}/calibration_and_imaging.py (100%) diff --git a/demos/README.md b/demos/README.md new file mode 100644 index 00000000..fa19ca00 --- /dev/null +++ b/demos/README.md @@ -0,0 +1,16 @@ +# How to run the demos + +## Basic imaging with automatic weighting + +- Download the [data](https://www.philipp-arras.de/assets/CYG-ALL-2052-2MHZ.ms.tar.gz) and unpack it. +- Change the path under `[data]` in `cygnusa.cfg` to the path where the data is located. +- Run +```sh +python3 imaging_with_automatic_weighting.py cygnusa.cfg +``` +or, if you have `mpi4py` installed: + +``` sh +mpirun -np <ntasks> python3 imaging_with_automatic_weighting.py cygnusa.cfg +``` +which should speed up the computation. The number of threads used *per task* can be set in the configuration file `cygnusa.cfg` in the section `[technical]/nthreads`. The number threads multiplied by the number of MPI tasks should not exceed the number CPU cores available on the machine. diff --git a/demo/cygnusa.cfg b/demos/cygnusa.cfg similarity index 100% rename from demo/cygnusa.cfg rename to demos/cygnusa.cfg diff --git a/demo/imaging_with_automatic_weighting.py b/demos/imaging_with_automatic_weighting.py similarity index 100% rename from demo/imaging_with_automatic_weighting.py rename to demos/imaging_with_automatic_weighting.py diff --git a/demo/basic_polarization.py b/demos/old/basic_polarization.py similarity index 100% rename from demo/basic_polarization.py rename to demos/old/basic_polarization.py diff --git a/demo/calibration_and_imaging.py b/demos/old/calibration_and_imaging.py similarity index 100% rename from demo/calibration_and_imaging.py rename to demos/old/calibration_and_imaging.py -- GitLab