From 8e327d70436071b234ab0fd383cb1b40152d4f22 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Wed, 3 Feb 2021 11:36:14 +0100
Subject: [PATCH 01/37] Implement mosaicing response
---
demo/mosaicing.py | 104 +++++++++++++++++++++++++++
resolve/__init__.py | 7 +-
resolve/likelihood.py | 19 +++--
resolve/mosaicing/__init__.py | 1 +
resolve/mosaicing/sky_slicer.py | 124 ++++++++++++++++++++++++++++++++
resolve/observation.py | 2 +
resolve/primary_beam.py | 14 ++++
resolve/response.py | 11 +++
8 files changed, 273 insertions(+), 9 deletions(-)
create mode 100644 demo/mosaicing.py
create mode 100644 resolve/mosaicing/__init__.py
create mode 100644 resolve/mosaicing/sky_slicer.py
diff --git a/demo/mosaicing.py b/demo/mosaicing.py
new file mode 100644
index 00000000..3d0ee31d
--- /dev/null
+++ b/demo/mosaicing.py
@@ -0,0 +1,104 @@
+# SPDX-License-Identifier: GPL-3.0-or-later
+# Copyright(C) 2019-2020 Max-Planck-Society
+# Author: Philipp Arras
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+import nifty7 as ift
+import resolve as rve
+
+
+def main():
+ rve.set_nthreads(1)
+ rve.set_wgridding(False)
+ diffusefluxlevel = 20
+
+ interobs = {
+ f"extended{ii}": rve.Observation.load(
+ f"/data/mosaicing/extended/extended_field{ii}.npz"
+ )
+ for ii in range(5, 135)
+ # for ii in [50, 58, 59, 68, 69, 70, 88]
+ # for ii in [68, 69, 70]
+ }
+ rve.set_epsilon(1 / 10 / max(o.max_snr() for o in interobs.values()))
+
+ # Plot individual dirty images
+ for kk, vv in interobs.items():
+ npix = 500
+ fov = 1 * rve.ARCMIN2RAD
+ domain = ift.RGSpace([npix, npix], [fov / npix, fov / npix])
+ R = rve.StokesIResponse(vv, domain)
+ ift.single_plot(R.adjoint(vv.vis * vv.weight), name=f"dirty{kk}.png")
+
+ xs, ys = [], []
+ for vv in interobs.values():
+ xs.append(vv.direction.phase_center[0])
+ ys.append(vv.direction.phase_center[1])
+ xs, ys = np.array(xs), np.array(ys)
+ # Convention: The phase center is the coordinate of the [nx/2, ny/2] pixel
+ global_phase_center_casa = -2.7180339078174175, -0.5210930907116116
+ global_phase_center = (max(xs) + min(xs)) / 2, (max(ys) + min(ys)) / 2
+ margin = 2 * rve.ARCMIN2RAD
+ plt.scatter(xs, ys, label="Extended")
+ plt.scatter(*global_phase_center_casa, label="Phase center CASA")
+ plt.scatter(*global_phase_center, label="Phase center resolve")
+ xfov = 2 * max(abs(xs - global_phase_center[0])) + 2 * margin
+ yfov = 2 * max(abs(ys - global_phase_center[1])) + 2 * margin
+ print(f"Field of view: {xfov/rve.ARCMIN2RAD:.2f}' x {yfov/rve.ARCMIN2RAD:.2f}'")
+ plt.xlim([global_phase_center[0] - xfov / 2, global_phase_center[0] + xfov / 2])
+ plt.ylim([global_phase_center[1] - yfov / 2, global_phase_center[1] + yfov / 2])
+ plt.legend()
+ plt.savefig("pointings.png")
+ plt.close()
+
+ # Assert all frequencies are the same
+ freq = list(interobs.values())[0].freq
+ for vv in interobs.values():
+ assert np.all(vv.freq == freq)
+ extended_blockage = np.array([10.7, 0.75]) # m
+ extended_blockage *= np.mean(freq) / rve.SPEEDOFLIGHT
+ beam_func = lambda x: rve.alma_beam_func(*extended_blockage, x)
+ beam_directions = {
+ kk: rve.BeamDirection(
+ obs.direction.phase_center[0] - global_phase_center[0],
+ obs.direction.phase_center[1] - global_phase_center[1],
+ beam_func,
+ 0.05,
+ )
+ for kk, obs in interobs.items()
+ }
+
+ npix = np.array([1000, 1000])
+ fov = np.array([xfov, yfov])
+ dom = ift.RGSpace(npix, fov / npix)
+ logsky = ift.SimpleCorrelatedField(
+ dom, diffusefluxlevel, (1, 0.1), (5, 1), (1.2, 0.4), (0.2, 0.2), (-2, 0.5)
+ )
+ skyslicer = rve.SkySlicer(logsky.target, beam_directions)
+ ift.extra.check_linear_operator(skyslicer)
+ ift.single_plot(
+ skyslicer.adjoint(ift.full(skyslicer.target, 1.0)),
+ name="debug_skyslicer_adjoint.png",
+ )
+ lh = rve.ImagingLikelihood(interobs, skyslicer @ logsky.exp())
+
+ plotter = rve.Plotter("png", "plots")
+ plotter.add("logsky", logsky)
+ plotter.add("power spectrum logsky", logsky.power_spectrum)
+ plotter.add_histogram(
+ "normalized residuals (original weights)", lh.normalized_residual
+ )
+ ham = ift.StandardHamiltonian(lh)
+ fld = 0.1 * ift.from_random(lh.domain)
+ state = rve.MinimizationState(fld, [])
+ mini = ift.NewtonCG(ift.GradientNormController(name="newton", iteration_limit=10))
+ state = rve.simple_minimize(ham, state.mean, 0, mini)
+ plotter.plot("stage1", state)
+ state.save("stage1")
+ # state = rve.MinimizationState.load("stage1")
+
+
+if __name__ == "__main__":
+ main()
diff --git a/resolve/__init__.py b/resolve/__init__.py
index fadc803f..5b1b670e 100644
--- a/resolve/__init__.py
+++ b/resolve/__init__.py
@@ -1,9 +1,12 @@
+from .antenna_positions import AntennaPositions
from .calibration import calibration_distribution
from .constants import *
+from .direction import *
from .fits import field2fits
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 *
from .ms_import import ms2observations, ms_n_spectral_windows
@@ -16,8 +19,8 @@ from .multi_frequency.operators import (
from .observation import Observation, tmin_tmax, unique_antennas, unique_times
from .plotter import MfPlotter, Plotter
from .points import PointInserter
-from .polarization import polarization_matrix_exponential
-from .primary_beam import vla_beam
+from .polarization import Polarization, polarization_matrix_exponential
+from .primary_beam import *
from .response import MfResponse, ResponseDistributor, StokesIResponse, SingleResponse
from .simple_operators import *
from .util import (
diff --git a/resolve/likelihood.py b/resolve/likelihood.py
index 44c3ad5e..e588b5e9 100644
--- a/resolve/likelihood.py
+++ b/resolve/likelihood.py
@@ -8,7 +8,7 @@ import nifty7 as ift
from .observation import Observation
from .response import FullPolResponse, MfResponse, StokesIResponse
-from .util import my_assert_isinstance, my_asserteq, my_assert
+from .util import my_assert, my_assert_isinstance, my_asserteq
def _get_mask(observation):
@@ -91,7 +91,7 @@ def ImagingLikelihood(
Parameters
----------
- observation : Observation
+ observation : Observation or dict(Observation)
Observation object from which observation.vis and potentially
observation.weight is used for computing the likelihood.
@@ -109,11 +109,12 @@ def ImagingLikelihood(
calibration_operator : Operator
Optional. Target needs to be the same as observation.vis.
"""
- my_assert_isinstance(observation, Observation)
my_assert_isinstance(sky_operator, ift.Operator)
sdom = sky_operator.target
- if isinstance(sdom, ift.MultiDomain):
+ 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."
@@ -121,14 +122,18 @@ def ImagingLikelihood(
else:
R = FullPolResponse(observation, sky_operator.target)
else:
- if len(sdom.shape) == 3:
+ 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
-
if inverse_covariance_operator is None:
- return _build_gauss_lh_nres(model_data, observation.vis, observation.weight)
+ if mosaicing:
+ vis = ift.MultiField.from_dict({kk: o.vis for kk, o in observation.items()})
+ weight = ift.MultiField.from_dict({kk: o.weight for kk, o in observation.items()})
+ else:
+ vis, weight = observation.vis, observation.weight
+ return _build_gauss_lh_nres(model_data, vis, weight)
return _varcov(observation, model_data, inverse_covariance_operator)
diff --git a/resolve/mosaicing/__init__.py b/resolve/mosaicing/__init__.py
new file mode 100644
index 00000000..6802616c
--- /dev/null
+++ b/resolve/mosaicing/__init__.py
@@ -0,0 +1 @@
+from .sky_slicer import *
diff --git a/resolve/mosaicing/sky_slicer.py b/resolve/mosaicing/sky_slicer.py
new file mode 100644
index 00000000..5fcf4253
--- /dev/null
+++ b/resolve/mosaicing/sky_slicer.py
@@ -0,0 +1,124 @@
+# SPDX-License-Identifier: GPL-3.0-or-later
+# Copyright(C) 2020 Max-Planck-Society
+# Author: Philipp Arras
+
+import numpy as np
+
+import nifty7 as ift
+
+
+class SkySlicer(ift.LinearOperator):
+ """Maps from the total sky domain to individual sky domains and applies the
+ primary beam pattern.
+
+ Parameters
+ ----------
+ domain : RGSpace
+ Two-dimensional RG-Space which serves as domain. The distances are
+ in pseudo-radian.
+
+ beam_directions : dict(key: BeamDirection)
+ Dictionary of BeamDirection that contains the information of the
+ directions of the different observations and the beam pattern.
+ """
+
+ def __init__(self, domain, beam_directions):
+ self._bd = dict(beam_directions)
+ self._domain = ift.makeDomain(domain)
+ t, b, s = {}, {}, {}
+ for kk, vv in self._bd.items():
+ t[kk], s[kk], b[kk] = vv.slice_target(self._domain)
+ self._beams = b
+ self._slices = s
+ self._target = ift.makeDomain(t)
+ 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 = {}
+ for kk in self._target.keys():
+ res[kk] = self._beams[kk].val * x[self._slices[kk]]
+ else:
+ res = np.zeros(self._domain.shape)
+ for kk, xx in x.items():
+ res[self._slices[kk]] += xx * self._beams[kk].val
+ return ift.makeField(self._tgt(mode), res)
+
+
+class BeamDirection:
+ """Represent direction information of one pointing.
+
+ Parameters
+ ----------
+ dx : float
+ Pointing offset in pseudo radian.
+ dy : float
+ Pointing offset in pseudo radian.
+ beam_func : function
+ Function that takes the two directions on the sky in pseudo radian and
+ returns beam pattern.
+ cutoff : float
+ Relative area of beam pattern that is cut off. 0 corresponds to no
+ cutoff at all. 1 corresponds to cut off everything.
+ """
+
+ def __init__(self, dx, dy, beam_func, cutoff):
+ self._dx, self._dy = float(dx), float(dy)
+ self._f = beam_func
+ self._cutoff = float(cutoff)
+ assert 0. <= cutoff < 1
+
+ def slice_target(self, domain):
+ """
+ Parameters
+ ----------
+ domain : RGSpace
+ Total sky domain.
+ """
+ dom = ift.makeDomain(domain)[0]
+ dst = np.array(dom.distances)
+ assert abs(self._dx) < dst[0] * dom.shape[0]
+ assert abs(self._dy) < dst[1] * dom.shape[1]
+ npix = np.array(domain.shape)
+
+ xs = np.linspace(0, max(npix * dst), num=4*max(npix))
+ # Technically not 100% correct since we integrate circles here.
+ # The actual cutoff is lower than the one that is induced by the
+ # following computation.
+ ys = self._f(xs)*xs
+ cond = np.cumsum(ys)/np.sum(ys)
+ np.testing.assert_allclose(cond[-1], 1.)
+ ind = np.searchsorted(cond, 1.-self._cutoff)
+ fov = 2*xs[ind]
+ patch_npix = fov/dst
+ patch_npix = (np.ceil(patch_npix/2)*2).astype(int)
+
+ shift = np.array([self._dx, self._dy])
+ patch_center_unrounded = np.array(npix) / 2 + shift / dst
+ ipix_patch_center = np.round(patch_center_unrounded).astype(int)
+ # Or maybe use ceil?
+ assert np.all(patch_npix % 2 == 0)
+ mi = (ipix_patch_center - patch_npix / 2).astype(int)
+ ma = mi + patch_npix
+
+ assert np.all(mi >= 0)
+ assert np.all(ma < npix)
+ slc = slice(mi[0], ma[0]), slice(mi[1], ma[1])
+ tgt = ift.RGSpace(patch_npix, dst)
+ assert tgt.shape[0] % 2 == 0
+ assert tgt.shape[1] % 2 == 0
+ test = np.empty(dom.shape)
+ assert test[slc].shape == tgt.shape
+
+ # Create coordinate field
+ mgrid = np.mgrid[: patch_npix[0], : patch_npix[1]].astype(float)
+ mgrid -= patch_npix[..., None, None] / 2
+ # FIXME Add subpixel offset
+ # subpixel_offset = ipix_patch_center - patch_center_unrounded
+ # assert np.all(subpixel_offset < 1.0)
+ # assert np.all(subpixel_offset >= 0.0)
+ mgrid *= dst[..., None, None]
+ beam = ift.makeField(tgt, self._f(np.linalg.norm(mgrid, axis=0)))
+ return tgt, slc, beam
diff --git a/resolve/observation.py b/resolve/observation.py
index 69735a55..6285292d 100644
--- a/resolve/observation.py
+++ b/resolve/observation.py
@@ -48,6 +48,8 @@ class Observation:
my_asserteq(vis.shape, (len(polarization), nrows, len(freq)))
my_asserteq(nrows, vis.shape[1])
my_assert(np.all(weight >= 0.0))
+ my_assert(np.all(np.isfinite(vis)))
+ my_assert(np.all(np.isfinite(weight)))
vis.flags.writeable = False
weight.flags.writeable = False
diff --git a/resolve/primary_beam.py b/resolve/primary_beam.py
index 71f8ebaa..b4a186dd 100644
--- a/resolve/primary_beam.py
+++ b/resolve/primary_beam.py
@@ -3,6 +3,7 @@
# Author: Philipp Arras
import numpy as np
+import scipy.special as sc
import nifty7 as ift
@@ -176,3 +177,16 @@ def vla_beam(domain, freq):
beam = rweight * upper + (1 - rweight) * lower
beam[beam < 0] = 0
return ift.makeOp(ift.makeField(dom, beam))
+
+
+def alma_beam_func(D, d, x):
+ assert isinstance(x, np.ndarray)
+ assert x.ndim < 3
+ assert np.max(np.abs(x)) < np.pi/np.sqrt(2)
+ b = d / D
+ x = np.pi * D * x
+ mask = x == 0.0
+ x[mask] = 1
+ sol = 2 / (x * (1 - b ** 2)) * (sc.jn(1, x) - b * sc.jn(1, x * b))
+ sol[mask] = 1
+ return sol**2
diff --git a/resolve/response.py b/resolve/response.py
index 89dbe374..f44931bc 100644
--- a/resolve/response.py
+++ b/resolve/response.py
@@ -2,6 +2,9 @@
# Copyright(C) 2019-2020 Max-Planck-Society
# Author: Philipp Arras
+from functools import reduce
+from operator import add
+
import numpy as np
from ducc0.wgridder.experimental import dirty2vis, vis2dirty
@@ -14,6 +17,14 @@ from .util import my_assert, my_assert_isinstance, my_asserteq
def StokesIResponse(observation, domain):
+ if isinstance(observation, dict):
+ # TODO Possibly add subpixel offset here
+ d = ift.MultiDomain.make(domain)
+ res = (
+ StokesIResponse(o, d[kk]).ducktape(kk).ducktape_left(kk)
+ for kk, o in observation.items()
+ )
+ return reduce(add, res)
my_assert_isinstance(observation, Observation)
domain = ift.DomainTuple.make(domain)
my_asserteq(len(domain), 1)
--
GitLab
From 4b04558ce77c67a1e2269163a266fbedaad34418 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Thu, 4 Feb 2021 17:02:33 +0100
Subject: [PATCH 02/37] Add alma data script
---
misc/convert_alma_data.py | 61 +++++++++++++++++++++++++++++++++++++++
1 file changed, 61 insertions(+)
create mode 100644 misc/convert_alma_data.py
diff --git a/misc/convert_alma_data.py b/misc/convert_alma_data.py
new file mode 100644
index 00000000..da327660
--- /dev/null
+++ b/misc/convert_alma_data.py
@@ -0,0 +1,61 @@
+import numpy as np
+import resolve as rve
+from os.path import join
+from casacore.tables import table
+
+f = "extended.ms"
+_CASACORE_TABLE_CFG = {"readonly": True, "ack": False}
+
+spectral_window = 0
+channel_slice = slice(1027, 1030)
+
+with table(join(f, "SPECTRAL_WINDOW"), **_CASACORE_TABLE_CFG) as t:
+ freq = np.array(t.getcol("CHAN_FREQ"))[spectral_window]
+freq = freq[channel_slice]
+
+with table(join(f, "FIELD"), **_CASACORE_TABLE_CFG) as t:
+ field_directions = np.squeeze(t.getcol("DELAY_DIR"))
+assert field_directions.ndim == 2
+assert field_directions.shape[1] == 2
+
+with table(f, **_CASACORE_TABLE_CFG) as t:
+ wgt = np.array(t.getcol("WEIGHT"))[:, None]
+ wgt = np.repeat(wgt, len(freq), axis=1)
+ flag = np.array(t.getcol("FLAG"))[:, channel_slice]
+ wgt[flag] = 0.
+ datadescid = np.array(t.getcol("DATA_DESC_ID"))
+ wgt[datadescid != spectral_window] = 0.
+ del(flag, spectral_window)
+
+ vis = np.array(t.getcol("DATA"))[:, channel_slice]
+ uvw = np.array(t.getcol("UVW"))
+ field = np.array(t.getcol("FIELD_ID"))
+
+ # Average polarization
+ assert vis.shape[2] == 2
+ assert wgt.shape[2] == 2
+ vis = np.sum(wgt*vis, axis=-1)
+ wgt = np.sum(wgt, axis=2)
+ vis /= wgt
+
+for ifield in set(field):
+ ww = wgt.copy()
+ ww[field != ifield] = 0.
+
+ # Clean up rows
+ ind = np.any(ww != 0., axis=1)
+ myuvw = uvw[ind]
+ mywgt = wgt[ind]
+ myvis = vis[ind]
+
+ # Clean up freqs
+ assert ww.ndim == 2
+ ind = np.any(mywgt != 0., axis=0)
+ mywgt = mywgt[:, ind]
+ myvis = myvis[:, ind]
+ myfreq = freq[ind]
+ assert len(myfreq) == 3
+ print(ifield, myuvw.shape, myvis.shape, mywgt.shape, myfreq.shape)
+
+ obs = rve.Observation(rve.AntennaPositions(myuvw), myvis[None], mywgt[None], rve.Polarization.trivial(), myfreq, rve.Direction(field_directions[ifield], 2000))
+ obs.save(f"extended_field{ifield}.npz", compress=False)
--
GitLab
From d3ee790e6bab7c36dc5efe2e7a1315d86491c571 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Thu, 4 Feb 2021 17:09:48 +0100
Subject: [PATCH 03/37] Save np arrays as contiguous
---
misc/convert_alma_data.py | 32 +++++++++++++++++++++-----------
1 file changed, 21 insertions(+), 11 deletions(-)
diff --git a/misc/convert_alma_data.py b/misc/convert_alma_data.py
index da327660..e24bc4b6 100644
--- a/misc/convert_alma_data.py
+++ b/misc/convert_alma_data.py
@@ -1,9 +1,12 @@
-import numpy as np
-import resolve as rve
+import sys
from os.path import join
+
+import numpy as np
from casacore.tables import table
-f = "extended.ms"
+import resolve as rve
+
+f = sys.argv[1]
_CASACORE_TABLE_CFG = {"readonly": True, "ack": False}
spectral_window = 0
@@ -22,10 +25,10 @@ with table(f, **_CASACORE_TABLE_CFG) as t:
wgt = np.array(t.getcol("WEIGHT"))[:, None]
wgt = np.repeat(wgt, len(freq), axis=1)
flag = np.array(t.getcol("FLAG"))[:, channel_slice]
- wgt[flag] = 0.
+ wgt[flag] = 0.0
datadescid = np.array(t.getcol("DATA_DESC_ID"))
- wgt[datadescid != spectral_window] = 0.
- del(flag, spectral_window)
+ wgt[datadescid != spectral_window] = 0.0
+ del (flag, spectral_window)
vis = np.array(t.getcol("DATA"))[:, channel_slice]
uvw = np.array(t.getcol("UVW"))
@@ -34,28 +37,35 @@ with table(f, **_CASACORE_TABLE_CFG) as t:
# Average polarization
assert vis.shape[2] == 2
assert wgt.shape[2] == 2
- vis = np.sum(wgt*vis, axis=-1)
+ vis = np.sum(wgt * vis, axis=-1)
wgt = np.sum(wgt, axis=2)
vis /= wgt
for ifield in set(field):
ww = wgt.copy()
- ww[field != ifield] = 0.
+ ww[field != ifield] = 0.0
# Clean up rows
- ind = np.any(ww != 0., axis=1)
+ ind = np.any(ww != 0.0, axis=1)
myuvw = uvw[ind]
mywgt = wgt[ind]
myvis = vis[ind]
# Clean up freqs
assert ww.ndim == 2
- ind = np.any(mywgt != 0., axis=0)
+ ind = np.any(mywgt != 0.0, axis=0)
mywgt = mywgt[:, ind]
myvis = myvis[:, ind]
myfreq = freq[ind]
assert len(myfreq) == 3
print(ifield, myuvw.shape, myvis.shape, mywgt.shape, myfreq.shape)
- obs = rve.Observation(rve.AntennaPositions(myuvw), myvis[None], mywgt[None], rve.Polarization.trivial(), myfreq, rve.Direction(field_directions[ifield], 2000))
+ obs = rve.Observation(
+ rve.AntennaPositions(np.ascontiguousarray(myuvw)),
+ np.ascontiguousarray(myvis[None]),
+ np.ascontiguousarray(mywgt[None]),
+ rve.Polarization.trivial(),
+ myfreq,
+ rve.Direction(field_directions[ifield], 2000),
+ )
obs.save(f"extended_field{ifield}.npz", compress=False)
--
GitLab
From 96ccf287aa5b25b9e1e7b5d89136f120c9f6ff09 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Fri, 5 Feb 2021 11:03:36 +0100
Subject: [PATCH 04/37] Add debugging info
---
demo/mosaicing.py | 18 ++++++++----------
resolve/response.py | 21 +++++++++++++++++++++
2 files changed, 29 insertions(+), 10 deletions(-)
diff --git a/demo/mosaicing.py b/demo/mosaicing.py
index 3d0ee31d..0e50be28 100644
--- a/demo/mosaicing.py
+++ b/demo/mosaicing.py
@@ -10,7 +10,7 @@ import resolve as rve
def main():
- rve.set_nthreads(1)
+ rve.set_nthreads(2)
rve.set_wgridding(False)
diffusefluxlevel = 20
@@ -24,14 +24,6 @@ def main():
}
rve.set_epsilon(1 / 10 / max(o.max_snr() for o in interobs.values()))
- # Plot individual dirty images
- for kk, vv in interobs.items():
- npix = 500
- fov = 1 * rve.ARCMIN2RAD
- domain = ift.RGSpace([npix, npix], [fov / npix, fov / npix])
- R = rve.StokesIResponse(vv, domain)
- ift.single_plot(R.adjoint(vv.vis * vv.weight), name=f"dirty{kk}.png")
-
xs, ys = [], []
for vv in interobs.values():
xs.append(vv.direction.phase_center[0])
@@ -57,6 +49,7 @@ def main():
freq = list(interobs.values())[0].freq
for vv in interobs.values():
assert np.all(vv.freq == freq)
+ print("Frequencies:", freq)
extended_blockage = np.array([10.7, 0.75]) # m
extended_blockage *= np.mean(freq) / rve.SPEEDOFLIGHT
beam_func = lambda x: rve.alma_beam_func(*extended_blockage, x)
@@ -70,7 +63,7 @@ def main():
for kk, obs in interobs.items()
}
- npix = np.array([1000, 1000])
+ npix = np.array([1500, 1500])
fov = np.array([xfov, yfov])
dom = ift.RGSpace(npix, fov / npix)
logsky = ift.SimpleCorrelatedField(
@@ -84,6 +77,11 @@ def main():
)
lh = rve.ImagingLikelihood(interobs, skyslicer @ logsky.exp())
+ R = rve.StokesIResponse(interobs, skyslicer.target) @ skyslicer
+ vis = ift.MultiField.from_dict({kk: o.vis for kk, o in interobs.items()})
+ weight = ift.MultiField.from_dict({kk: o.weight for kk, o in interobs.items()})
+ ift.single_plot(R.adjoint(vis*weight), name="total_dirty.png")
+
plotter = rve.Plotter("png", "plots")
plotter.add("logsky", logsky)
plotter.add("power spectrum logsky", logsky.power_spectrum)
diff --git a/resolve/response.py b/resolve/response.py
index f44931bc..43fade8b 100644
--- a/resolve/response.py
+++ b/resolve/response.py
@@ -10,6 +10,7 @@ from ducc0.wgridder.experimental import dirty2vis, vis2dirty
import nifty7 as ift
+from .constants import SPEEDOFLIGHT
from .global_config import epsilon, nthreads, wgridding
from .multi_frequency.irg_space import IRGSpace
from .observation import Observation
@@ -238,6 +239,7 @@ class SingleResponse(ift.LinearOperator):
self._target_dtype = np.complex64 if single_precision else np.complex128
self._domain_dtype = np.float32 if single_precision else np.float64
self._verbt, self._verbadj = verbose, verbose
+ self._ofac = None
def apply(self, x, mode):
self._check_input(x, mode)
@@ -250,6 +252,9 @@ class SingleResponse(ift.LinearOperator):
args1 = {"dirty": x}
if self._verbt:
args1["verbosity"] = True
+ print(
+ f"\nINFO: Oversampling factors in response: {self.oversampling_factors()}\n"
+ )
self._verbt = False
f = dirty2vis
# FIXME Use vis_out keyword of wgridder
@@ -263,6 +268,9 @@ class SingleResponse(ift.LinearOperator):
}
if self._verbadj:
args1["verbosity"] = True
+ print(
+ f"\nINFO: Oversampling factors in response: {self.oversampling_factors()}\n"
+ )
self._verbadj = False
f = vis2dirty
res = ift.makeField(self._tgt(mode), f(**self._args, **args1) * self._vol)
@@ -270,3 +278,16 @@ class SingleResponse(ift.LinearOperator):
res.dtype, self._target_dtype if mode == self.TIMES else self._domain_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
--
GitLab
From 3a6d869c8a29ddb41e71576a2afe0c22128abcc2 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Fri, 5 Feb 2021 11:06:08 +0100
Subject: [PATCH 05/37] Improve data import
---
misc/convert_alma_data.py | 5 +++--
1 file changed, 3 insertions(+), 2 deletions(-)
diff --git a/misc/convert_alma_data.py b/misc/convert_alma_data.py
index e24bc4b6..643c60e1 100644
--- a/misc/convert_alma_data.py
+++ b/misc/convert_alma_data.py
@@ -1,5 +1,5 @@
import sys
-from os.path import join
+from os.path import join, split, splitext
import numpy as np
from casacore.tables import table
@@ -7,6 +7,7 @@ from casacore.tables import table
import resolve as rve
f = sys.argv[1]
+outfbase = split(splitext(f)[0])[1]
_CASACORE_TABLE_CFG = {"readonly": True, "ack": False}
spectral_window = 0
@@ -68,4 +69,4 @@ for ifield in set(field):
myfreq,
rve.Direction(field_directions[ifield], 2000),
)
- obs.save(f"extended_field{ifield}.npz", compress=False)
+ obs.save(f"{outfbase}_field{ifield}.npz", compress=False)
--
GitLab
From 72d135130b7cf6cd5c22eb270424e24827e0ebe8 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Fri, 5 Feb 2021 14:05:31 +0100
Subject: [PATCH 06/37] Add crazy offset in data import
---
misc/convert_alma_data.py | 12 +++++++++++-
1 file changed, 11 insertions(+), 1 deletion(-)
diff --git a/misc/convert_alma_data.py b/misc/convert_alma_data.py
index 643c60e1..41649cda 100644
--- a/misc/convert_alma_data.py
+++ b/misc/convert_alma_data.py
@@ -11,10 +11,20 @@ outfbase = split(splitext(f)[0])[1]
_CASACORE_TABLE_CFG = {"readonly": True, "ack": False}
spectral_window = 0
-channel_slice = slice(1027, 1030)
with table(join(f, "SPECTRAL_WINDOW"), **_CASACORE_TABLE_CFG) as t:
freq = np.array(t.getcol("CHAN_FREQ"))[spectral_window]
+
+minfreq = 230172900000.
+if outfbase == "compact":
+ minfreq -= 27146364.99822998
+maxfreq = minfreq * 1.000015
+a = freq[::-1]
+assert np.all(np.diff(a) >= 0)
+ind1 = freq.size - np.searchsorted(a, minfreq) + 1
+ind0 = freq.size - np.searchsorted(a, maxfreq) + 1
+channel_slice = slice(ind0, ind1)
+print(channel_slice)
freq = freq[channel_slice]
with table(join(f, "FIELD"), **_CASACORE_TABLE_CFG) as t:
--
GitLab
From 770b86465cdcbeb4ad2275ac53c9bd6709316341 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Fri, 5 Feb 2021 14:07:10 +0100
Subject: [PATCH 07/37] Fixup
---
misc/convert_alma_data.py | 1 -
1 file changed, 1 deletion(-)
diff --git a/misc/convert_alma_data.py b/misc/convert_alma_data.py
index 41649cda..2002f089 100644
--- a/misc/convert_alma_data.py
+++ b/misc/convert_alma_data.py
@@ -68,7 +68,6 @@ for ifield in set(field):
mywgt = mywgt[:, ind]
myvis = myvis[:, ind]
myfreq = freq[ind]
- assert len(myfreq) == 3
print(ifield, myuvw.shape, myvis.shape, mywgt.shape, myfreq.shape)
obs = rve.Observation(
--
GitLab
From 422e4adf99342a8d8c1c80a9bbaa77a52a13c66e Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Fri, 5 Feb 2021 14:08:16 +0100
Subject: [PATCH 08/37] Implement mosaicing response 9/n
---
demo/mosaicing.py | 35 +++++++++++++++------------------
resolve/mosaicing/sky_slicer.py | 2 +-
2 files changed, 17 insertions(+), 20 deletions(-)
diff --git a/demo/mosaicing.py b/demo/mosaicing.py
index 0e50be28..9f9e801e 100644
--- a/demo/mosaicing.py
+++ b/demo/mosaicing.py
@@ -14,14 +14,18 @@ def main():
rve.set_wgridding(False)
diffusefluxlevel = 20
- interobs = {
- f"extended{ii}": rve.Observation.load(
- f"/data/mosaicing/extended/extended_field{ii}.npz"
- )
+ interobs_extended = {
+ f"extended{ii}": rve.Observation.load(f"/data/mosaicing/extended_field{ii}.npz")
for ii in range(5, 135)
# for ii in [50, 58, 59, 68, 69, 70, 88]
# for ii in [68, 69, 70]
}
+ interobs_compact = {
+ f"compact{ii}": rve.Observation.load(f"/data/mosaicing/compact_field{ii}.npz")
+ for ii in range(5, 153)
+ }
+ interobs = {**interobs_compact, **interobs_extended}
+
rve.set_epsilon(1 / 10 / max(o.max_snr() for o in interobs.values()))
xs, ys = [], []
@@ -45,25 +49,18 @@ def main():
plt.savefig("pointings.png")
plt.close()
- # Assert all frequencies are the same
- freq = list(interobs.values())[0].freq
- for vv in interobs.values():
- assert np.all(vv.freq == freq)
- print("Frequencies:", freq)
- extended_blockage = np.array([10.7, 0.75]) # m
- extended_blockage *= np.mean(freq) / rve.SPEEDOFLIGHT
- beam_func = lambda x: rve.alma_beam_func(*extended_blockage, x)
- beam_directions = {
- kk: rve.BeamDirection(
+ beam_directions = {}
+ for kk, obs in interobs.items():
+ extended_blockage = np.array([10.7, 0.75]) # m
+ extended_blockage *= np.mean(obs.freq) / rve.SPEEDOFLIGHT
+ beam_directions[kk] = rve.BeamDirection(
obs.direction.phase_center[0] - global_phase_center[0],
obs.direction.phase_center[1] - global_phase_center[1],
- beam_func,
+ lambda x: rve.alma_beam_func(*extended_blockage, x),
0.05,
)
- for kk, obs in interobs.items()
- }
- npix = np.array([1500, 1500])
+ npix = np.array([1800, 1800])
fov = np.array([xfov, yfov])
dom = ift.RGSpace(npix, fov / npix)
logsky = ift.SimpleCorrelatedField(
@@ -80,7 +77,7 @@ def main():
R = rve.StokesIResponse(interobs, skyslicer.target) @ skyslicer
vis = ift.MultiField.from_dict({kk: o.vis for kk, o in interobs.items()})
weight = ift.MultiField.from_dict({kk: o.weight for kk, o in interobs.items()})
- ift.single_plot(R.adjoint(vis*weight), name="total_dirty.png")
+ plt.imsave("total_dirty.png", R.adjoint(vis * weight).val.T, origin="lower")
plotter = rve.Plotter("png", "plots")
plotter.add("logsky", logsky)
diff --git a/resolve/mosaicing/sky_slicer.py b/resolve/mosaicing/sky_slicer.py
index 5fcf4253..7eb75e4d 100644
--- a/resolve/mosaicing/sky_slicer.py
+++ b/resolve/mosaicing/sky_slicer.py
@@ -39,7 +39,7 @@ class SkySlicer(ift.LinearOperator):
if mode == self.TIMES:
res = {}
for kk in self._target.keys():
- res[kk] = self._beams[kk].val * x[self._slices[kk]]
+ res[kk] = x[self._slices[kk]] * self._beams[kk].val
else:
res = np.zeros(self._domain.shape)
for kk, xx in x.items():
--
GitLab
From febfe0bebea66a04371908d96a4b7d94a7d816c7 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Fri, 5 Feb 2021 14:08:40 +0100
Subject: [PATCH 09/37] Crucial fix in pointing conventions
---
resolve/mosaicing/sky_slicer.py | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/resolve/mosaicing/sky_slicer.py b/resolve/mosaicing/sky_slicer.py
index 7eb75e4d..3f377890 100644
--- a/resolve/mosaicing/sky_slicer.py
+++ b/resolve/mosaicing/sky_slicer.py
@@ -95,7 +95,7 @@ class BeamDirection:
patch_npix = fov/dst
patch_npix = (np.ceil(patch_npix/2)*2).astype(int)
- shift = np.array([self._dx, self._dy])
+ shift = np.array([-self._dx, self._dy])
patch_center_unrounded = np.array(npix) / 2 + shift / dst
ipix_patch_center = np.round(patch_center_unrounded).astype(int)
# Or maybe use ceil?
--
GitLab
From 1e0f47033bb4972a45d84aa656c64bd811dca405 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Fri, 5 Feb 2021 14:55:50 +0100
Subject: [PATCH 10/37] Remove dirty hack
---
misc/convert_alma_data.py | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/misc/convert_alma_data.py b/misc/convert_alma_data.py
index 2002f089..7a0e0a6b 100644
--- a/misc/convert_alma_data.py
+++ b/misc/convert_alma_data.py
@@ -16,8 +16,8 @@ with table(join(f, "SPECTRAL_WINDOW"), **_CASACORE_TABLE_CFG) as t:
freq = np.array(t.getcol("CHAN_FREQ"))[spectral_window]
minfreq = 230172900000.
-if outfbase == "compact":
- minfreq -= 27146364.99822998
+#if outfbase == "compact":
+ #minfreq -= 27146364.99822998
maxfreq = minfreq * 1.000015
a = freq[::-1]
assert np.all(np.diff(a) >= 0)
--
GitLab
From 8c4586c334f85d77f9029c375e15bb94342c9003 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Fri, 5 Feb 2021 16:16:09 +0100
Subject: [PATCH 11/37] Set up first working MGVI reconstruction
---
demo/mosaicing.py | 44 +++++++++++++++++++++++++++-----------------
1 file changed, 27 insertions(+), 17 deletions(-)
diff --git a/demo/mosaicing.py b/demo/mosaicing.py
index 9f9e801e..1cdba5e8 100644
--- a/demo/mosaicing.py
+++ b/demo/mosaicing.py
@@ -9,6 +9,12 @@ import nifty7 as ift
import resolve as rve
+@rve.onlymaster
+def imsave(fname, field):
+ plt.imsave(fname, field.val.T, origin="lower")
+ plt.close()
+
+
def main():
rve.set_nthreads(2)
rve.set_wgridding(False)
@@ -36,7 +42,7 @@ def main():
# Convention: The phase center is the coordinate of the [nx/2, ny/2] pixel
global_phase_center_casa = -2.7180339078174175, -0.5210930907116116
global_phase_center = (max(xs) + min(xs)) / 2, (max(ys) + min(ys)) / 2
- margin = 2 * rve.ARCMIN2RAD
+ margin = 1 * rve.ARCMIN2RAD
plt.scatter(xs, ys, label="Extended")
plt.scatter(*global_phase_center_casa, label="Phase center CASA")
plt.scatter(*global_phase_center, label="Phase center resolve")
@@ -45,8 +51,11 @@ def main():
print(f"Field of view: {xfov/rve.ARCMIN2RAD:.2f}' x {yfov/rve.ARCMIN2RAD:.2f}'")
plt.xlim([global_phase_center[0] - xfov / 2, global_phase_center[0] + xfov / 2])
plt.ylim([global_phase_center[1] - yfov / 2, global_phase_center[1] + yfov / 2])
+ plt.gca().invert_xaxis()
+ plt.gca().set_aspect("equal")
plt.legend()
- plt.savefig("pointings.png")
+ if rve.mpi.master:
+ plt.savefig("pointings.png")
plt.close()
beam_directions = {}
@@ -57,41 +66,42 @@ def main():
obs.direction.phase_center[0] - global_phase_center[0],
obs.direction.phase_center[1] - global_phase_center[1],
lambda x: rve.alma_beam_func(*extended_blockage, x),
- 0.05,
+ 0.1,
)
npix = np.array([1800, 1800])
fov = np.array([xfov, yfov])
dom = ift.RGSpace(npix, fov / npix)
logsky = ift.SimpleCorrelatedField(
- dom, diffusefluxlevel, (1, 0.1), (5, 1), (1.2, 0.4), (0.2, 0.2), (-2, 0.5)
+ dom, diffusefluxlevel, (1, 0.1), (1.5, 1), (1.2, 0.4), (0.2, 0.2), (-1.5, 0.5)
)
skyslicer = rve.SkySlicer(logsky.target, beam_directions)
ift.extra.check_linear_operator(skyslicer)
- ift.single_plot(
- skyslicer.adjoint(ift.full(skyslicer.target, 1.0)),
- name="debug_skyslicer_adjoint.png",
- )
lh = rve.ImagingLikelihood(interobs, skyslicer @ logsky.exp())
R = rve.StokesIResponse(interobs, skyslicer.target) @ skyslicer
vis = ift.MultiField.from_dict({kk: o.vis for kk, o in interobs.items()})
weight = ift.MultiField.from_dict({kk: o.weight for kk, o in interobs.items()})
- plt.imsave("total_dirty.png", R.adjoint(vis * weight).val.T, origin="lower")
+
+ imsave("total_dirty.png", R.adjoint(vis * weight))
+ imsave("skyslicer_adjoint.png", skyslicer.adjoint(ift.full(skyslicer.target, 1.0)))
plotter = rve.Plotter("png", "plots")
plotter.add("logsky", logsky)
plotter.add("power spectrum logsky", logsky.power_spectrum)
- plotter.add_histogram(
- "normalized residuals (original weights)", lh.normalized_residual
- )
- ham = ift.StandardHamiltonian(lh)
+ # plotter.add_histogram(
+ # "normalized residuals (original weights)", lh.normalized_residual
+ # )
+ ham = ift.StandardHamiltonian(lh, ift.AbsDeltaEnergyController(0.5, 3, 300))
fld = 0.1 * ift.from_random(lh.domain)
state = rve.MinimizationState(fld, [])
- mini = ift.NewtonCG(ift.GradientNormController(name="newton", iteration_limit=10))
- state = rve.simple_minimize(ham, state.mean, 0, mini)
- plotter.plot("stage1", state)
- state.save("stage1")
+ mini = ift.NewtonCG(ift.GradientNormController(name="newton", iteration_limit=5))
+ # TODO Add minisanity to resolve
+ # TODO Fix histogram plots
+ for ii in range(20):
+ state = rve.simple_minimize(ham, state.mean, 3, mini)
+ plotter.plot(f"iter{ii}", state)
+ # state.save("stage1")
# state = rve.MinimizationState.load("stage1")
--
GitLab
From ae87a91d3e9cb056aadc35bf7e17ff0ce0f941f5 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Sun, 7 Feb 2021 12:47:37 +0100
Subject: [PATCH 12/37] Minor
---
demo/mosaicing.py | 6 +++---
1 file changed, 3 insertions(+), 3 deletions(-)
diff --git a/demo/mosaicing.py b/demo/mosaicing.py
index 1cdba5e8..1cd74980 100644
--- a/demo/mosaicing.py
+++ b/demo/mosaicing.py
@@ -88,11 +88,12 @@ def main():
plotter = rve.Plotter("png", "plots")
plotter.add("logsky", logsky)
+ plotter.add("sky", logsky.exp(), cmap="inferno_r")
plotter.add("power spectrum logsky", logsky.power_spectrum)
# plotter.add_histogram(
# "normalized residuals (original weights)", lh.normalized_residual
# )
- ham = ift.StandardHamiltonian(lh, ift.AbsDeltaEnergyController(0.5, 3, 300))
+ ham = ift.StandardHamiltonian(lh, ift.AbsDeltaEnergyController(0.5, 3, 300, name="Sampling"))
fld = 0.1 * ift.from_random(lh.domain)
state = rve.MinimizationState(fld, [])
mini = ift.NewtonCG(ift.GradientNormController(name="newton", iteration_limit=5))
@@ -101,8 +102,7 @@ def main():
for ii in range(20):
state = rve.simple_minimize(ham, state.mean, 3, mini)
plotter.plot(f"iter{ii}", state)
- # state.save("stage1")
- # state = rve.MinimizationState.load("stage1")
+ state.save(f"iter{ii}")
if __name__ == "__main__":
--
GitLab
From 564d8e7046ff320265eedb5e94ba4324b90569d2 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Sun, 7 Feb 2021 17:55:39 +0100
Subject: [PATCH 13/37] Start integrating single dish data
---
misc/convert_alma_data.py | 175 ++++++++++++++++++++++--------------
resolve/observation.py | 181 +++++++++++++++++++++++++++-----------
2 files changed, 237 insertions(+), 119 deletions(-)
diff --git a/misc/convert_alma_data.py b/misc/convert_alma_data.py
index 7a0e0a6b..65d7e1c5 100644
--- a/misc/convert_alma_data.py
+++ b/misc/convert_alma_data.py
@@ -1,81 +1,124 @@
-import sys
+# SPDX-License-Identifier: GPL-3.0-or-later
+# Copyright(C) 2021 Max-Planck-Society
+# Author: Philipp Arras
+
from os.path import join, split, splitext
import numpy as np
-from casacore.tables import table
+from casacore.tables import table, taql
import resolve as rve
-f = sys.argv[1]
-outfbase = split(splitext(f)[0])[1]
-_CASACORE_TABLE_CFG = {"readonly": True, "ack": False}
-
-spectral_window = 0
-
-with table(join(f, "SPECTRAL_WINDOW"), **_CASACORE_TABLE_CFG) as t:
- freq = np.array(t.getcol("CHAN_FREQ"))[spectral_window]
-
-minfreq = 230172900000.
-#if outfbase == "compact":
- #minfreq -= 27146364.99822998
-maxfreq = minfreq * 1.000015
-a = freq[::-1]
-assert np.all(np.diff(a) >= 0)
-ind1 = freq.size - np.searchsorted(a, minfreq) + 1
-ind0 = freq.size - np.searchsorted(a, maxfreq) + 1
-channel_slice = slice(ind0, ind1)
-print(channel_slice)
-freq = freq[channel_slice]
-
-with table(join(f, "FIELD"), **_CASACORE_TABLE_CFG) as t:
- field_directions = np.squeeze(t.getcol("DELAY_DIR"))
-assert field_directions.ndim == 2
-assert field_directions.shape[1] == 2
-
-with table(f, **_CASACORE_TABLE_CFG) as t:
- wgt = np.array(t.getcol("WEIGHT"))[:, None]
- wgt = np.repeat(wgt, len(freq), axis=1)
- flag = np.array(t.getcol("FLAG"))[:, channel_slice]
- wgt[flag] = 0.0
- datadescid = np.array(t.getcol("DATA_DESC_ID"))
- wgt[datadescid != spectral_window] = 0.0
- del (flag, spectral_window)
-
- vis = np.array(t.getcol("DATA"))[:, channel_slice]
- uvw = np.array(t.getcol("UVW"))
- field = np.array(t.getcol("FIELD_ID"))
-
- # Average polarization
- assert vis.shape[2] == 2
- assert wgt.shape[2] == 2
- vis = np.sum(wgt * vis, axis=-1)
- wgt = np.sum(wgt, axis=2)
- vis /= wgt
-
-for ifield in set(field):
- ww = wgt.copy()
- ww[field != ifield] = 0.0
+def cleanup(weight, vis, freq, uvw=None):
+ assert weight.ndim == 2
+ assert vis.shape == weight.shape
+ assert freq.ndim==1
+ assert freq.size == vis.shape[1]
+ assert uvw.shape[0] == vis.shape[0]
# Clean up rows
- ind = np.any(ww != 0.0, axis=1)
- myuvw = uvw[ind]
- mywgt = wgt[ind]
+ ind = np.any(weight != 0.0, axis=1)
+ mywgt = weight[ind]
myvis = vis[ind]
+ if uvw is not None:
+ myuvw = uvw[ind]
# Clean up freqs
- assert ww.ndim == 2
ind = np.any(mywgt != 0.0, axis=0)
mywgt = mywgt[:, ind]
myvis = myvis[:, ind]
myfreq = freq[ind]
- print(ifield, myuvw.shape, myvis.shape, mywgt.shape, myfreq.shape)
-
- obs = rve.Observation(
- rve.AntennaPositions(np.ascontiguousarray(myuvw)),
- np.ascontiguousarray(myvis[None]),
- np.ascontiguousarray(mywgt[None]),
- rve.Polarization.trivial(),
- myfreq,
- rve.Direction(field_directions[ifield], 2000),
- )
- obs.save(f"{outfbase}_field{ifield}.npz", compress=False)
+ if uvw is None:
+ return mywgt, myvis, myfreq
+ return mywgt, myvis, myfreq, myuvw
+
+
+_CASACORE_TABLE_CFG = {"readonly": True, "ack": False}
+minfreq = 230172900000.0
+# if outfbase == "compact":
+# minfreq -= 27146364.99822998
+maxfreq = minfreq * 1.000015
+
+for f, spectral_window, single_dish, fieldid in [
+ # ("extended.ms", 0, False, None),
+ # ("compact.ms", 0, False, None),
+ ("totalpower.ms", 17, True, 1)
+]:
+ outfbase = split(splitext(f)[0])[1]
+
+ with table(join(f, "SPECTRAL_WINDOW"), **_CASACORE_TABLE_CFG) as t:
+ freq = taql("select CHAN_FREQ from $t where ROWID()=$spectral_window").getcol(
+ "CHAN_FREQ"
+ )
+ freq = np.squeeze(freq)
+
+ if not np.all(np.diff(freq) >= 0):
+ a = freq[::-1]
+ assert np.all(np.diff(a) >= 0)
+ ind1 = freq.size - np.searchsorted(a, minfreq) + 1
+ ind0 = freq.size - np.searchsorted(a, maxfreq) + 1
+ else:
+ assert np.all(np.diff(freq) >= 0)
+ ind0 = np.searchsorted(freq, minfreq)
+ ind1 = np.searchsorted(freq, maxfreq)
+ channel_slice = slice(ind0, ind1)
+ freq = freq[channel_slice]
+
+ if not single_dish:
+ with table(join(f, "FIELD"), **_CASACORE_TABLE_CFG) as t:
+ field_directions = np.squeeze(t.getcol("DELAY_DIR"))
+ assert field_directions.ndim == 2
+ assert field_directions.shape[1] == 2
+
+ with table(f, **_CASACORE_TABLE_CFG) as t:
+ sel = lambda x: np.array(
+ taql(f"select {x} from $t where DATA_DESC_ID=$spectral_window").getcol(x)
+ )
+ wgt = sel("WEIGHT")[:, None]
+ wgt = np.repeat(wgt, len(freq), axis=1)
+ flag = sel("FLAG")[:, channel_slice]
+ wgt[flag] = 0.0
+ del flag
+
+ vis = sel("FLOAT_DATA" if single_dish else "DATA")[:, channel_slice]
+ # Average polarization
+ assert vis.shape[2] == 2
+ assert wgt.shape[2] == 2
+ vis = np.sum(wgt * vis, axis=-1)
+ wgt = np.sum(wgt, axis=2)
+ vis /= wgt
+ field = sel("FIELD_ID")
+ if single_dish:
+ wgt[field != fieldid] = 0.0
+ del field
+ else:
+ uvw = sel("UVW")
+
+ if single_dish:
+ with table(f, readonly=True) as t:
+ pointings = taql("select DIRECTION from $t::POINTING where ROWID() in [select from $t where DATA_DESC_ID=$spectral_window giving [ROWID()]]").getcol("DIRECTION")
+ mywgt, myvis, myfreq, mypointings = cleanup(wgt, vis, freq, pointings)
+ obs = rve.SingleDishObservation(
+ np.ascontiguousarray(mypointings),
+ np.ascontiguousarray(myvis[None]),
+ np.ascontiguousarray(mywgt[None]),
+ rve.Polarization.trivial(),
+ myfreq,
+ )
+ obs.save(f"{outfbase}.npz", compress=False)
+ else:
+ for ifield in set(field):
+ ww = wgt.copy()
+ ww[field != ifield] = 0.0
+ mywgt, myvis, myfreq, myuvw = cleanup(ww, vis, freq, uvw)
+ print(ifield, myuvw.shape, myvis.shape, mywgt.shape, myfreq.shape)
+
+ obs = rve.Observation(
+ rve.AntennaPositions(np.ascontiguousarray(myuvw)),
+ np.ascontiguousarray(myvis[None]),
+ np.ascontiguousarray(mywgt[None]),
+ rve.Polarization.trivial(),
+ myfreq,
+ rve.Direction(field_directions[ifield], 2000),
+ )
+ obs.save(f"{outfbase}_field{ifield}.npz", compress=False)
diff --git a/resolve/observation.py b/resolve/observation.py
index 6285292d..c89c1253 100644
--- a/resolve/observation.py
+++ b/resolve/observation.py
@@ -1,5 +1,5 @@
# SPDX-License-Identifier: GPL-3.0-or-later
-# Copyright(C) 2019-2020 Max-Planck-Society
+# Copyright(C) 2019-2021 Max-Planck-Society
# Author: Philipp Arras
import numpy as np
@@ -14,7 +14,133 @@ from .polarization import Polarization
from .util import compare_attributes, my_assert, my_assert_isinstance, my_asserteq
-class Observation:
+class _Observation:
+ @property
+ def vis(self):
+ dom = [ift.UnstructuredDomain(ss) for ss in self._vis.shape]
+ return ift.makeField(dom, self._vis)
+
+ @property
+ def weight(self):
+ dom = [ift.UnstructuredDomain(ss) for ss in self._weight.shape]
+ return ift.makeField(dom, self._weight)
+
+ @property
+ def freq(self):
+ return self._freq
+
+ @property
+ def polarization(self):
+ return self._polarization
+
+ @property
+ def direction(self):
+ return self._direction
+
+ @property
+ def npol(self):
+ return self._vis.shape[0]
+
+ @property
+ def nrow(self):
+ return self._vis.shape[1]
+
+ @property
+ def nfreq(self):
+ return self._vis.shape[2]
+
+ def apply_flags(self, arr):
+ return arr[self._weight != 0.0]
+
+ @property
+ def flags(self):
+ return self._weight == 0.0
+
+ @property
+ def mask(self):
+ return self._weight > 0.0
+
+ def max_snr(self):
+ return np.max(np.abs(self.apply_flags(self._vis * np.sqrt(self._weight))))
+
+ def fraction_useful(self):
+ return self.apply_flags(self._weight).size / self._weight.size
+
+
+class SingleDishObservation(_Observation):
+ def __init__(self):
+ self._weight
+ self._vis
+ self._direction
+ self._polarization
+ self._freq
+
+ @onlymaster
+ def save(self, file_name, compress):
+ raise NotImplementedError
+ # dct = dict(
+ # vis=self._vis,
+ # 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:
+ # vv = np.array([])
+ # dct[f"antpos{ii}"] = vv
+ # f = np.savez_compressed if compress else np.savez
+ # f(file_name, **dct)
+
+ @staticmethod
+ def load(file_name):
+ raise NotImplementedError
+ # dct = dict(np.load(file_name))
+ # antpos = []
+ # for ii in range(4):
+ # val = dct[f"antpos{ii}"]
+ # if val.size == 0:
+ # val = None
+ # antpos.append(val)
+ # pol = Polarization.from_list(dct["polarization"])
+ # direction = Direction.from_list(dct["direction"])
+ # return Observation(
+ # AntennaPositions.from_list(antpos),
+ # dct["vis"],
+ # dct["weight"],
+ # pol,
+ # dct["freq"],
+ # direction,
+ # )
+
+ def __eq__(self, other):
+ raise NotImplementedError
+ # if not isinstance(other, Observation):
+ # return False
+ # if (
+ # self._vis.dtype != other._vis.dtype
+ # or self._weight.dtype != other._weight.dtype
+ # ):
+ # return False
+ # return compare_attributes(
+ # self,
+ # other,
+ # ("_direction", "_polarization", "_freq", "_antpos", "_vis", "_weight"),
+ # )
+
+ def __getitem__(self, slc):
+ raise NotImplementedError
+ # return Observation(
+ # self._antpos[slc],
+ # self._vis[:, slc],
+ # self._weight[:, slc],
+ # self._polarization,
+ # self._freq,
+ # self._direction,
+ # )
+
+
+class Observation(_Observation):
"""Observation data
This class contains all the data and information about an observation.
@@ -61,23 +187,6 @@ class Observation:
self._freq = freq
self._direction = direction
- def apply_flags(self, arr):
- return arr[self._weight != 0.0]
-
- @property
- def flags(self):
- return self._weight == 0.0
-
- @property
- def mask(self):
- return self._weight > 0.0
-
- def max_snr(self):
- return np.max(np.abs(self.apply_flags(self._vis * np.sqrt(self._weight))))
-
- def fraction_useful(self):
- return self.apply_flags(self._weight).size / self._weight.size
-
@onlymaster
def save(self, file_name, compress):
dct = dict(
@@ -172,40 +281,6 @@ class Observation:
uvlen = np.linalg.norm(self.uvw, axis=1)
return np.outer(uvlen, self._freq / SPEEDOFLIGHT)
- @property
- def vis(self):
- dom = [ift.UnstructuredDomain(ss) for ss in self._vis.shape]
- return ift.makeField(dom, self._vis)
-
- @property
- def weight(self):
- dom = [ift.UnstructuredDomain(ss) for ss in self._weight.shape]
- return ift.makeField(dom, self._weight)
-
- @property
- def freq(self):
- return self._freq
-
- @property
- def polarization(self):
- return self._polarization
-
- @property
- def direction(self):
- return self._direction
-
- @property
- def npol(self):
- return self._vis.shape[0]
-
- @property
- def nrow(self):
- return self._vis.shape[1]
-
- @property
- def nfreq(self):
- return self._vis.shape[2]
-
def tmin_tmax(*args):
"""
--
GitLab
From beea25d0300b2aedaba0f3c50a5af405aeb16845 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Sun, 7 Feb 2021 21:50:41 +0100
Subject: [PATCH 14/37] Add SingleDishObservation
---
misc/convert_alma_data.py | 13 ++++-
resolve/__init__.py | 2 +-
resolve/direction.py | 39 ++++++++++++-
resolve/observation.py | 117 ++++++++++++++++++--------------------
4 files changed, 105 insertions(+), 66 deletions(-)
diff --git a/misc/convert_alma_data.py b/misc/convert_alma_data.py
index 65d7e1c5..0962f982 100644
--- a/misc/convert_alma_data.py
+++ b/misc/convert_alma_data.py
@@ -13,7 +13,7 @@ import resolve as rve
def cleanup(weight, vis, freq, uvw=None):
assert weight.ndim == 2
assert vis.shape == weight.shape
- assert freq.ndim==1
+ assert freq.ndim == 1
assert freq.size == vis.shape[1]
assert uvw.shape[0] == vis.shape[0]
# Clean up rows
@@ -96,10 +96,15 @@ for f, spectral_window, single_dish, fieldid in [
if single_dish:
with table(f, readonly=True) as t:
- pointings = taql("select DIRECTION from $t::POINTING where ROWID() in [select from $t where DATA_DESC_ID=$spectral_window giving [ROWID()]]").getcol("DIRECTION")
+ pointings = taql(
+ "select DIRECTION from $t::POINTING where ROWID() in [select from $t where DATA_DESC_ID=$spectral_window giving [ROWID()]]"
+ ).getcol("DIRECTION")
mywgt, myvis, myfreq, mypointings = cleanup(wgt, vis, freq, pointings)
+ assert mypointings.shape[1] == 1
+ mypointings = mypointings[:, 0]
+ mypointings = rve.Directions(mypointings, 2000)
obs = rve.SingleDishObservation(
- np.ascontiguousarray(mypointings),
+ mypointings,
np.ascontiguousarray(myvis[None]),
np.ascontiguousarray(mywgt[None]),
rve.Polarization.trivial(),
@@ -122,3 +127,5 @@ for f, spectral_window, single_dish, fieldid in [
rve.Direction(field_directions[ifield], 2000),
)
obs.save(f"{outfbase}_field{ifield}.npz", compress=False)
+ check = rve.Observation.load(f"{outfbase}_field{ifield}.npz")
+ assert check == obs
diff --git a/resolve/__init__.py b/resolve/__init__.py
index 5b1b670e..b74b7fda 100644
--- a/resolve/__init__.py
+++ b/resolve/__init__.py
@@ -16,7 +16,7 @@ from .multi_frequency.operators import (
MfWeightingInterpolation,
WienerIntegrations,
)
-from .observation import Observation, tmin_tmax, unique_antennas, unique_times
+from .observation import *
from .plotter import MfPlotter, Plotter
from .points import PointInserter
from .polarization import Polarization, polarization_matrix_exponential
diff --git a/resolve/direction.py b/resolve/direction.py
index cc148054..a6b7b2fc 100644
--- a/resolve/direction.py
+++ b/resolve/direction.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 .util import compare_attributes, my_asserteq
@@ -42,3 +42,40 @@ class Direction:
if not isinstance(other, Direction):
return False
return compare_attributes(self, other, ("_pc", "_e"))
+
+
+class Directions:
+ def __init__(self, phase_centers, equinox):
+ assert phase_centers.ndim == 2
+ assert phase_centers.shape[1] == 2
+ self._pc = phase_centers
+ self._e = float(equinox)
+
+ @property
+ def phase_centers(self):
+ return self._pc
+
+ @property
+ def equinox(self):
+ return self._e
+
+ def __repr__(self):
+ return f"Directions({self._pc}, equinox={self._e})"
+
+ def to_list(self):
+ return [self._pc, self._e]
+
+ def __len__(self):
+ return self._pc.shape[0]
+
+ @staticmethod
+ def from_list(lst):
+ return Direction(lst[0], lst[1])
+
+ def __eq__(self, other):
+ if not isinstance(other, Direction):
+ return False
+ return compare_attributes(self, other, ("_pc", "_e"))
+
+ def __getitem__(self, slc):
+ return Directions(self._pc[slc], self._e)
diff --git a/resolve/observation.py b/resolve/observation.py
index c89c1253..44355333 100644
--- a/resolve/observation.py
+++ b/resolve/observation.py
@@ -8,7 +8,7 @@ import nifty7 as ift
from .antenna_positions import AntennaPositions
from .constants import SPEEDOFLIGHT
-from .direction import Direction
+from .direction import Direction, Directions
from .mpi import onlymaster
from .polarization import Polarization
from .util import compare_attributes, my_assert, my_assert_isinstance, my_asserteq
@@ -68,76 +68,71 @@ class _Observation:
class SingleDishObservation(_Observation):
- def __init__(self):
- self._weight
- self._vis
- self._direction
- self._polarization
- self._freq
+ def __init__(self, pointings, data, weight, polarization, freq):
+ my_assert_isinstance(pointings, Directions)
+ my_assert_isinstance(polarization, Polarization)
+ my_assert(data.dtype in [np.float32, np.float64])
+ nrows = len(pointings)
+ my_asserteq(weight.shape, data.shape)
+ my_asserteq(data.shape, (len(polarization), nrows, len(freq)))
+ my_asserteq(nrows, data.shape[1])
+
+ data.flags.writeable = False
+ weight.flags.writeable = False
+
+ my_assert(np.all(weight >= 0.0))
+ my_assert(np.all(np.isfinite(data)))
+ my_assert(np.all(np.isfinite(weight)))
+
+ self._pointings = pointings
+ self._vis = data
+ self._weight = weight
+ self._polarization = polarization
+ self._freq = freq
@onlymaster
def save(self, file_name, compress):
- raise NotImplementedError
- # dct = dict(
- # vis=self._vis,
- # 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:
- # vv = np.array([])
- # dct[f"antpos{ii}"] = vv
- # f = np.savez_compressed if compress else np.savez
- # f(file_name, **dct)
+ p = self._pointings.to_list()
+ dct = dict(
+ vis=self._vis,
+ weight=self._weight,
+ freq=self._freq,
+ polarization=self._polarization.to_list(),
+ pointings0=p[0],
+ pointings1=p[1],
+ )
+ f = np.savez_compressed if compress else np.savez
+ f(file_name, **dct)
@staticmethod
def load(file_name):
- raise NotImplementedError
- # dct = dict(np.load(file_name))
- # antpos = []
- # for ii in range(4):
- # val = dct[f"antpos{ii}"]
- # if val.size == 0:
- # val = None
- # antpos.append(val)
- # pol = Polarization.from_list(dct["polarization"])
- # direction = Direction.from_list(dct["direction"])
- # return Observation(
- # AntennaPositions.from_list(antpos),
- # dct["vis"],
- # dct["weight"],
- # pol,
- # dct["freq"],
- # direction,
- # )
+ dct = dict(np.load(file_name))
+ pol = Polarization.from_list(dct["polarization"])
+ pointings = Directions.from_list(dct["pointings0"], dct["pointings1"])
+ return SingleDishObservation(
+ pointings, dct["vis"], dct["weight"], pol, dct["freq"]
+ )
def __eq__(self, other):
- raise NotImplementedError
- # if not isinstance(other, Observation):
- # return False
- # if (
- # self._vis.dtype != other._vis.dtype
- # or self._weight.dtype != other._weight.dtype
- # ):
- # return False
- # return compare_attributes(
- # self,
- # other,
- # ("_direction", "_polarization", "_freq", "_antpos", "_vis", "_weight"),
- # )
+ if not isinstance(other, Observation):
+ return False
+ if (
+ self._vis.dtype != other._vis.dtype
+ or self._weight.dtype != other._weight.dtype
+ ):
+ return False
+ return compare_attributes(
+ self, other, ("_polarization", "_freq", "_pointings", "_vis", "_weight")
+ )
def __getitem__(self, slc):
- raise NotImplementedError
- # return Observation(
- # self._antpos[slc],
- # self._vis[:, slc],
- # self._weight[:, slc],
- # self._polarization,
- # self._freq,
- # self._direction,
- # )
+ return SingleDishObservation(
+ self._pointings[slc],
+ self._vis[:, slc],
+ self._weight[:, slc],
+ self._polarization,
+ self._freq,
+ )
class Observation(_Observation):
--
GitLab
From f295d5bc6504710953c1271c04f7df89d2aeac16 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Sun, 7 Feb 2021 22:04:00 +0100
Subject: [PATCH 15/37] Load single dish data
---
demo/mosaicing.py | 22 +++++++++++++++++++---
resolve/direction.py | 2 +-
resolve/observation.py | 6 +++++-
3 files changed, 25 insertions(+), 5 deletions(-)
diff --git a/demo/mosaicing.py b/demo/mosaicing.py
index 1cd74980..803472fe 100644
--- a/demo/mosaicing.py
+++ b/demo/mosaicing.py
@@ -31,6 +31,7 @@ def main():
for ii in range(5, 153)
}
interobs = {**interobs_compact, **interobs_extended}
+ sdobs = {"sd0": rve.SingleDishObservation.load("/data/mosaicing/totalpower.npz")}
rve.set_epsilon(1 / 10 / max(o.max_snr() for o in interobs.values()))
@@ -38,17 +39,32 @@ def main():
for vv in interobs.values():
xs.append(vv.direction.phase_center[0])
ys.append(vv.direction.phase_center[1])
+ for vv in sdobs.values():
+ foo, bar = vv.pointings.phase_centers.T
+ xs.extend(foo)
+ ys.extend(bar)
xs, ys = np.array(xs), np.array(ys)
+
# Convention: The phase center is the coordinate of the [nx/2, ny/2] pixel
global_phase_center_casa = -2.7180339078174175, -0.5210930907116116
global_phase_center = (max(xs) + min(xs)) / 2, (max(ys) + min(ys)) / 2
margin = 1 * rve.ARCMIN2RAD
- plt.scatter(xs, ys, label="Extended")
- plt.scatter(*global_phase_center_casa, label="Phase center CASA")
- plt.scatter(*global_phase_center, label="Phase center resolve")
xfov = 2 * max(abs(xs - global_phase_center[0])) + 2 * margin
yfov = 2 * max(abs(ys - global_phase_center[1])) + 2 * margin
print(f"Field of view: {xfov/rve.ARCMIN2RAD:.2f}' x {yfov/rve.ARCMIN2RAD:.2f}'")
+
+ xs, ys = [], []
+ for vv in interobs.values():
+ xs.append(vv.direction.phase_center[0])
+ ys.append(vv.direction.phase_center[1])
+ xs, ys = np.array(xs), np.array(ys)
+ plt.scatter(xs, ys, label="Interferometer")
+ for kk, vv in sdobs.items():
+ plt.scatter(
+ *vv.pointings.phase_centers.T, label="Single dish " + kk, alpha=0.2, s=1
+ )
+ plt.scatter(*global_phase_center_casa, label="Phase center CASA")
+ plt.scatter(*global_phase_center, label="Phase center resolve")
plt.xlim([global_phase_center[0] - xfov / 2, global_phase_center[0] + xfov / 2])
plt.ylim([global_phase_center[1] - yfov / 2, global_phase_center[1] + yfov / 2])
plt.gca().invert_xaxis()
diff --git a/resolve/direction.py b/resolve/direction.py
index a6b7b2fc..5d2acdad 100644
--- a/resolve/direction.py
+++ b/resolve/direction.py
@@ -70,7 +70,7 @@ class Directions:
@staticmethod
def from_list(lst):
- return Direction(lst[0], lst[1])
+ return Directions(lst[0], lst[1])
def __eq__(self, other):
if not isinstance(other, Direction):
diff --git a/resolve/observation.py b/resolve/observation.py
index 44355333..3d2a0a28 100644
--- a/resolve/observation.py
+++ b/resolve/observation.py
@@ -108,7 +108,7 @@ class SingleDishObservation(_Observation):
def load(file_name):
dct = dict(np.load(file_name))
pol = Polarization.from_list(dct["polarization"])
- pointings = Directions.from_list(dct["pointings0"], dct["pointings1"])
+ pointings = Directions.from_list([dct["pointings0"], dct["pointings1"]])
return SingleDishObservation(
pointings, dct["vis"], dct["weight"], pol, dct["freq"]
)
@@ -134,6 +134,10 @@ class SingleDishObservation(_Observation):
self._freq,
)
+ @property
+ def pointings(self):
+ return self._pointings
+
class Observation(_Observation):
"""Observation data
--
GitLab
From 607171f44da67f3b2abfeecb50eecb93c3e825f7 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Mon, 8 Feb 2021 10:13:47 +0100
Subject: [PATCH 16/37] Tweak beam
---
demo/mosaicing.py | 5 ++---
resolve/primary_beam.py | 12 +++++++-----
2 files changed, 9 insertions(+), 8 deletions(-)
diff --git a/demo/mosaicing.py b/demo/mosaicing.py
index 803472fe..adffda54 100644
--- a/demo/mosaicing.py
+++ b/demo/mosaicing.py
@@ -76,12 +76,11 @@ def main():
beam_directions = {}
for kk, obs in interobs.items():
- extended_blockage = np.array([10.7, 0.75]) # m
- extended_blockage *= np.mean(obs.freq) / rve.SPEEDOFLIGHT
+ assert np.min(obs.freq)/np.max(obs.freq) > 0.99
beam_directions[kk] = rve.BeamDirection(
obs.direction.phase_center[0] - global_phase_center[0],
obs.direction.phase_center[1] - global_phase_center[1],
- lambda x: rve.alma_beam_func(*extended_blockage, x),
+ lambda x: rve.alma_beam_func(10.7, 0.75, np.mean(obs.freq), x),
0.1,
)
diff --git a/resolve/primary_beam.py b/resolve/primary_beam.py
index b4a186dd..f627f5e2 100644
--- a/resolve/primary_beam.py
+++ b/resolve/primary_beam.py
@@ -7,7 +7,7 @@ import scipy.special as sc
import nifty7 as ift
-from .constants import ARCMIN2RAD
+from .constants import ARCMIN2RAD, SPEEDOFLIGHT
from .util import my_assert
@@ -179,14 +179,16 @@ def vla_beam(domain, freq):
return ift.makeOp(ift.makeField(dom, beam))
-def alma_beam_func(D, d, x):
+def alma_beam_func(D, d, freq, x):
assert isinstance(x, np.ndarray)
assert x.ndim < 3
- assert np.max(np.abs(x)) < np.pi/np.sqrt(2)
+ assert np.max(np.abs(x)) < np.pi / np.sqrt(2)
+
+ a = freq / SPEEDOFLIGHT
b = d / D
- x = np.pi * D * x
+ x = np.pi * a * D * x
mask = x == 0.0
x[mask] = 1
sol = 2 / (x * (1 - b ** 2)) * (sc.jn(1, x) - b * sc.jn(1, x * b))
sol[mask] = 1
- return sol**2
+ return sol ** 2
--
GitLab
From d9a9411b0867e883c36bc479c0df9ce6127a268e Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Mon, 8 Feb 2021 11:54:45 +0100
Subject: [PATCH 17/37] Implement single dish likelihood
---
demo/mosaicing.py | 63 ++++++++++++++++++++++++++-------
resolve/likelihood.py | 30 ++++++++++++++++
resolve/mosaicing/__init__.py | 1 +
resolve/mosaicing/sky_slicer.py | 1 +
4 files changed, 82 insertions(+), 13 deletions(-)
diff --git a/demo/mosaicing.py b/demo/mosaicing.py
index adffda54..85327700 100644
--- a/demo/mosaicing.py
+++ b/demo/mosaicing.py
@@ -1,8 +1,10 @@
# SPDX-License-Identifier: GPL-3.0-or-later
-# Copyright(C) 2019-2020 Max-Planck-Society
+# Copyright(C) 2019-2021 Max-Planck-Society
# Author: Philipp Arras
import numpy as np
+from functools import reduce
+from operator import add
import matplotlib.pyplot as plt
import nifty7 as ift
@@ -32,9 +34,12 @@ def main():
}
interobs = {**interobs_compact, **interobs_extended}
sdobs = {"sd0": rve.SingleDishObservation.load("/data/mosaicing/totalpower.npz")}
+ assert len(set(sdobs.keys()) & set(interobs.keys())) == 0
rve.set_epsilon(1 / 10 / max(o.max_snr() for o in interobs.values()))
+ # Compute phase center and define domain
+ # Convention: The phase center is the coordinate of the [nx/2, ny/2] pixel
xs, ys = [], []
for vv in interobs.values():
xs.append(vv.direction.phase_center[0])
@@ -44,15 +49,18 @@ def main():
xs.extend(foo)
ys.extend(bar)
xs, ys = np.array(xs), np.array(ys)
-
- # Convention: The phase center is the coordinate of the [nx/2, ny/2] pixel
global_phase_center_casa = -2.7180339078174175, -0.5210930907116116
global_phase_center = (max(xs) + min(xs)) / 2, (max(ys) + min(ys)) / 2
margin = 1 * rve.ARCMIN2RAD
xfov = 2 * max(abs(xs - global_phase_center[0])) + 2 * margin
yfov = 2 * max(abs(ys - global_phase_center[1])) + 2 * margin
print(f"Field of view: {xfov/rve.ARCMIN2RAD:.2f}' x {yfov/rve.ARCMIN2RAD:.2f}'")
+ npix = np.array([1800, 1800])
+ fov = np.array([xfov, yfov])
+ dom = ift.RGSpace(npix, fov / npix)
+ # End compute phase center and define domain
+ # Plot pointings
xs, ys = [], []
for vv in interobs.values():
xs.append(vv.direction.phase_center[0])
@@ -73,33 +81,62 @@ def main():
if rve.mpi.master:
plt.savefig("pointings.png")
plt.close()
+ # End plot pointings
+
+ # Single dish likelihood
+ single_dish_response = []
+ for kk, oo in sdobs.items():
+ assert np.min(oo.freq) / np.max(oo.freq) > 0.99
+ R = rve.SingleDishResponse(
+ oo,
+ dom,
+ lambda x: rve.alma_beam_func(10.7, 0.75, np.mean(oo.freq), x),
+ global_phase_center,
+ )
+ single_dish_response.append(R.ducktape_left(kk))
+ single_dish_response = reduce(add, single_dish_response)
+ dsd = ift.MultiField.from_dict({kk: o.vis for kk, o in sdobs.items()})
+ invcovsd = ift.makeOp(
+ ift.MultiField.from_dict({kk: o.weight for kk, o in sdobs.items()})
+ )
+ lhsd = ift.GaussianEnergy(dsd, invcovsd) @ single_dish_response
+ imsave("sd_dirty.png", single_dish_response.adjoint(invcovsd(dsd)))
+ # End single dish likelihood
+ logsky = ift.SimpleCorrelatedField(
+ dom, diffusefluxlevel, (1, 0.1), (1.5, 1), (1.2, 0.4), (0.2, 0.2), (-1.5, 0.5)
+ )
+ sky = logsky.exp()
+ # p = ift.Plot()
+ # for _ in range(9):
+ # p.add(logsky(ift.from_random(logsky.domain)))
+ # p.output(name="prior.png")
+
+ # Interferometer likelihood
beam_directions = {}
for kk, obs in interobs.items():
- assert np.min(obs.freq)/np.max(obs.freq) > 0.99
+ assert np.min(obs.freq) / np.max(obs.freq) > 0.99
beam_directions[kk] = rve.BeamDirection(
obs.direction.phase_center[0] - global_phase_center[0],
obs.direction.phase_center[1] - global_phase_center[1],
lambda x: rve.alma_beam_func(10.7, 0.75, np.mean(obs.freq), x),
0.1,
)
-
- npix = np.array([1800, 1800])
- fov = np.array([xfov, yfov])
- dom = ift.RGSpace(npix, fov / npix)
- logsky = ift.SimpleCorrelatedField(
- dom, diffusefluxlevel, (1, 0.1), (1.5, 1), (1.2, 0.4), (0.2, 0.2), (-1.5, 0.5)
- )
skyslicer = rve.SkySlicer(logsky.target, beam_directions)
+ print(skyslicer.target)
ift.extra.check_linear_operator(skyslicer)
- lh = rve.ImagingLikelihood(interobs, skyslicer @ logsky.exp())
+ lhinter = rve.ImagingLikelihood(interobs, skyslicer)
+ # End interferometer likelihood
+ lh = (lhinter + lhsd) @ sky
+
+ # Plots
R = rve.StokesIResponse(interobs, skyslicer.target) @ skyslicer
vis = ift.MultiField.from_dict({kk: o.vis for kk, o in interobs.items()})
weight = ift.MultiField.from_dict({kk: o.weight for kk, o in interobs.items()})
-
imsave("total_dirty.png", R.adjoint(vis * weight))
imsave("skyslicer_adjoint.png", skyslicer.adjoint(ift.full(skyslicer.target, 1.0)))
+ # End plots
plotter = rve.Plotter("png", "plots")
plotter.add("logsky", logsky)
diff --git a/resolve/likelihood.py b/resolve/likelihood.py
index e588b5e9..c714f35c 100644
--- a/resolve/likelihood.py
+++ b/resolve/likelihood.py
@@ -63,6 +63,36 @@ def _varcov(observation, Rs, inverse_covariance_operator):
return _build_varcov_gauss_lh_nres(residual, inverse_covariance_operator, dtype)
+def SingleDishLikelihood(
+ observation,
+ sky_operator,
+ inverse_covariance_operator=None,
+ calibration_operator=None,
+):
+ if inverse_covariance_operator is not None or calibration_operator is not None:
+ raise NotImplementedError
+ mosaicing = isinstance(observation, dict)
+ if mosaicing:
+ vis = ift.MultiField.from_dict({kk: o.vis for kk, o in observation.items()})
+ weight = ift.MultiField.from_dict(
+ {kk: o.weight for kk, o in observation.items()}
+ )
+ single_dish_response = []
+ for kk, oo in observation.items():
+ assert np.min(oo.freq) / np.max(oo.freq) > 0.99
+ R = rve.SingleDishResponse(
+ oo,
+ sky_operator.target,
+ lambda x: rve.alma_beam_func(10.7, 0.75, np.mean(oo.freq), x),
+ global_phase_center,
+ )
+ single_dish_response.append(R.ducktape_left(kk))
+ from functools import reduce
+ from operator import add
+
+ single_dish_response = reduce(add, single_dish_response)
+
+
def ImagingLikelihood(
observation,
sky_operator,
diff --git a/resolve/mosaicing/__init__.py b/resolve/mosaicing/__init__.py
index 6802616c..a346d91e 100644
--- a/resolve/mosaicing/__init__.py
+++ b/resolve/mosaicing/__init__.py
@@ -1 +1,2 @@
from .sky_slicer import *
+from .single_dish_response import *
diff --git a/resolve/mosaicing/sky_slicer.py b/resolve/mosaicing/sky_slicer.py
index 3f377890..8ae42236 100644
--- a/resolve/mosaicing/sky_slicer.py
+++ b/resolve/mosaicing/sky_slicer.py
@@ -95,6 +95,7 @@ class BeamDirection:
patch_npix = fov/dst
patch_npix = (np.ceil(patch_npix/2)*2).astype(int)
+ # Convention: pointing convention (see also SingleDishResponse)
shift = np.array([-self._dx, self._dy])
patch_center_unrounded = np.array(npix) / 2 + shift / dst
ipix_patch_center = np.round(patch_center_unrounded).astype(int)
--
GitLab
From 421b4a0364a632b0a95508c792d737948a7b22e8 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Mon, 8 Feb 2021 12:08:48 +0100
Subject: [PATCH 18/37] Cosmetics
---
demo/mosaicing.py | 14 +++++---------
resolve/mosaicing/sky_slicer.py | 2 ++
2 files changed, 7 insertions(+), 9 deletions(-)
diff --git a/demo/mosaicing.py b/demo/mosaicing.py
index 85327700..1ddceabf 100644
--- a/demo/mosaicing.py
+++ b/demo/mosaicing.py
@@ -25,8 +25,6 @@ def main():
interobs_extended = {
f"extended{ii}": rve.Observation.load(f"/data/mosaicing/extended_field{ii}.npz")
for ii in range(5, 135)
- # for ii in [50, 58, 59, 68, 69, 70, 88]
- # for ii in [68, 69, 70]
}
interobs_compact = {
f"compact{ii}": rve.Observation.load(f"/data/mosaicing/compact_field{ii}.npz")
@@ -54,7 +52,8 @@ def main():
margin = 1 * rve.ARCMIN2RAD
xfov = 2 * max(abs(xs - global_phase_center[0])) + 2 * margin
yfov = 2 * max(abs(ys - global_phase_center[1])) + 2 * margin
- print(f"Field of view: {xfov/rve.ARCMIN2RAD:.2f}' x {yfov/rve.ARCMIN2RAD:.2f}'")
+ if rve.mpi.master:
+ print(f"Field of view: {xfov/rve.ARCMIN2RAD:.2f}' x {yfov/rve.ARCMIN2RAD:.2f}'")
npix = np.array([1800, 1800])
fov = np.array([xfov, yfov])
dom = ift.RGSpace(npix, fov / npix)
@@ -107,10 +106,6 @@ def main():
dom, diffusefluxlevel, (1, 0.1), (1.5, 1), (1.2, 0.4), (0.2, 0.2), (-1.5, 0.5)
)
sky = logsky.exp()
- # p = ift.Plot()
- # for _ in range(9):
- # p.add(logsky(ift.from_random(logsky.domain)))
- # p.output(name="prior.png")
# Interferometer likelihood
beam_directions = {}
@@ -123,7 +118,6 @@ def main():
0.1,
)
skyslicer = rve.SkySlicer(logsky.target, beam_directions)
- print(skyslicer.target)
ift.extra.check_linear_operator(skyslicer)
lhinter = rve.ImagingLikelihood(interobs, skyslicer)
# End interferometer likelihood
@@ -145,7 +139,9 @@ def main():
# plotter.add_histogram(
# "normalized residuals (original weights)", lh.normalized_residual
# )
- ham = ift.StandardHamiltonian(lh, ift.AbsDeltaEnergyController(0.5, 3, 300, name="Sampling"))
+ ham = ift.StandardHamiltonian(
+ lh, ift.AbsDeltaEnergyController(0.5, 5, 300, name="Sampling")
+ )
fld = 0.1 * ift.from_random(lh.domain)
state = rve.MinimizationState(fld, [])
mini = ift.NewtonCG(ift.GradientNormController(name="newton", iteration_limit=5))
diff --git a/resolve/mosaicing/sky_slicer.py b/resolve/mosaicing/sky_slicer.py
index 8ae42236..fb951280 100644
--- a/resolve/mosaicing/sky_slicer.py
+++ b/resolve/mosaicing/sky_slicer.py
@@ -27,7 +27,9 @@ class SkySlicer(ift.LinearOperator):
self._domain = ift.makeDomain(domain)
t, b, s = {}, {}, {}
for kk, vv in self._bd.items():
+ print("\r" + kk, end="")
t[kk], s[kk], b[kk] = vv.slice_target(self._domain)
+ print()
self._beams = b
self._slices = s
self._target = ift.makeDomain(t)
--
GitLab
From 6b3978a28f7fbb049f2b2bd78956d699ea2bab74 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Mon, 8 Feb 2021 12:08:56 +0100
Subject: [PATCH 19/37] Add file
---
resolve/mosaicing/single_dish_response.py | 32 +++++++++++++++++++++++
1 file changed, 32 insertions(+)
create mode 100644 resolve/mosaicing/single_dish_response.py
diff --git a/resolve/mosaicing/single_dish_response.py b/resolve/mosaicing/single_dish_response.py
new file mode 100644
index 00000000..d4aa2b87
--- /dev/null
+++ b/resolve/mosaicing/single_dish_response.py
@@ -0,0 +1,32 @@
+# SPDX-License-Identifier: GPL-3.0-or-later
+# Copyright(C) 2021 Max-Planck-Society
+# Author: Philipp Arras
+
+import numpy as np
+
+import nifty7 as ift
+
+from ..observation import SingleDishObservation
+
+
+def SingleDishResponse(observation, domain, beam_function, global_phase_center):
+ assert isinstance(observation, SingleDishObservation)
+ domain = ift.makeDomain(domain)
+ assert len(domain) == 1
+ codomain = domain[0].get_default_codomain()
+ kernel = codomain.get_conv_kernel_from_func(beam_function)
+ HT = ift.HartleyOperator(domain, codomain)
+ conv = HT.inverse @ ift.makeOp(kernel) @ HT.scale(domain.total_volume())
+ # FIXME Move into tests
+ fld = ift.from_random(conv.domain)
+ ift.extra.assert_allclose(conv(fld).integrate(), fld.integrate())
+ pc = observation.pointings.phase_centers.T - np.array(global_phase_center)[:, None]
+ pc = pc + (np.array(domain.shape)*np.array(domain[0].distances)/2)[:, None]
+ # Convention: pointing convention (see also BeamDirection)
+ pc[0] *= -1
+ interp = ift.LinearInterpolator(domain, pc)
+ bc = ift.ContractionOperator(observation.vis.domain, (0, 2)).adjoint
+ # NOTE The volume factor above `domain.total_volume()` and the volume factor
+ # below `domain[0].scalar_dvol` cancel each other. They are left in the
+ # code such that the convolution leaves the integral invariant.
+ return bc @ interp @ conv.scale(domain[0].scalar_dvol)
--
GitLab
From 71ce09c2dbf910a9ec05c95dfe6a524637a5fdec Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Mon, 8 Feb 2021 20:15:06 +0100
Subject: [PATCH 20/37] Restructure and add plots
---
demo/cfg_mosaicing.py | 99 ++++++++++++++++++++++++++++++++++
demo/mosaicing.py | 117 +----------------------------------------
demo/post_mosaicing.py | 112 +++++++++++++++++++++++++++++++++++++++
3 files changed, 213 insertions(+), 115 deletions(-)
create mode 100644 demo/cfg_mosaicing.py
create mode 100644 demo/post_mosaicing.py
diff --git a/demo/cfg_mosaicing.py b/demo/cfg_mosaicing.py
new file mode 100644
index 00000000..d34490d4
--- /dev/null
+++ b/demo/cfg_mosaicing.py
@@ -0,0 +1,99 @@
+# SPDX-License-Identifier: GPL-3.0-or-later
+# Copyright(C) 2019-2021 Max-Planck-Society
+# Author: Philipp Arras
+
+import numpy as np
+from functools import reduce
+from operator import add
+
+import nifty7 as ift
+import resolve as rve
+
+
+diffusefluxlevel = 20
+npix = np.array([1800, 1800])
+
+rve.set_nthreads(2)
+rve.set_wgridding(False)
+interobs_extended = {
+ f"extended{ii}": rve.Observation.load(f"/data/mosaicing/extended_field{ii}.npz")
+ for ii in range(5, 135)
+}
+interobs_compact = {
+ f"compact{ii}": rve.Observation.load(f"/data/mosaicing/compact_field{ii}.npz")
+ for ii in range(5, 153)
+}
+interobs = {**interobs_compact, **interobs_extended}
+sdobs = {"sd0": rve.SingleDishObservation.load("/data/mosaicing/totalpower.npz")}
+assert len(set(sdobs.keys()) & set(interobs.keys())) == 0
+rve.set_epsilon(1 / 10 / max(o.max_snr() for o in interobs.values()))
+# Compute phase center and define domain
+# Convention: The phase center is the coordinate of the [nx/2, ny/2] pixel
+xs, ys = [], []
+for vv in interobs.values():
+ xs.append(vv.direction.phase_center[0])
+ ys.append(vv.direction.phase_center[1])
+for vv in sdobs.values():
+ foo, bar = vv.pointings.phase_centers.T
+ xs.extend(foo)
+ ys.extend(bar)
+xs, ys = np.array(xs), np.array(ys)
+global_phase_center_casa = -2.7180339078174175, -0.5210930907116116
+global_phase_center = (max(xs) + min(xs)) / 2, (max(ys) + min(ys)) / 2
+margin = 1 * rve.ARCMIN2RAD
+xfov = 2 * max(abs(xs - global_phase_center[0])) + 2 * margin
+yfov = 2 * max(abs(ys - global_phase_center[1])) + 2 * margin
+fov = np.array([xfov, yfov])
+dom = ift.RGSpace(npix, fov / npix)
+# End compute phase center and define domain
+
+# Single dish likelihood
+responsesd = []
+for kk, oo in sdobs.items():
+ # Single dish response does not support multifrequency
+ assert np.min(oo.freq) / np.max(oo.freq) > 0.99
+ R = rve.SingleDishResponse(
+ oo,
+ dom,
+ lambda x: rve.alma_beam_func(10.7, 0.75, np.mean(oo.freq), x),
+ global_phase_center,
+ )
+ responsesd.append(R.ducktape_left(kk))
+responsesd = reduce(add, responsesd)
+dsd = ift.MultiField.from_dict({kk: o.vis for kk, o in sdobs.items()})
+invcovsd = ift.makeOp(
+ ift.MultiField.from_dict({kk: o.weight for kk, o in sdobs.items()})
+)
+lhsd = ift.GaussianEnergy(dsd, invcovsd) @ responsesd
+# End single dish likelihood
+
+# Sky model
+logsky = ift.SimpleCorrelatedField(
+ dom, diffusefluxlevel, (1, 0.1), (1.5, 1), (1.2, 0.4), (0.2, 0.2), (-1.5, 0.5)
+)
+sky = logsky.exp()
+# End sky model
+
+
+# Interferometer likelihood
+def get_beamdirections(obs):
+ beam_directions = {}
+ for kk, oo in obs.items():
+ assert np.min(oo.freq) / np.max(oo.freq) > 0.99
+ beam_directions[kk] = rve.BeamDirection(
+ oo.direction.phase_center[0] - global_phase_center[0],
+ oo.direction.phase_center[1] - global_phase_center[1],
+ lambda x: rve.alma_beam_func(10.7, 0.75, np.mean(oo.freq), x),
+ 0.1,
+ )
+ return beam_directions
+
+
+beamsc = get_beamdirections(interobs_compact)
+beamsext = get_beamdirections(interobs_extended)
+skyslicer_compact = rve.SkySlicer(logsky.target, beamsc)
+skyslicer_extended = rve.SkySlicer(logsky.target, beamsext)
+skyslicer = rve.SkySlicer(logsky.target, {**beamsc, **beamsext})
+ift.extra.check_linear_operator(skyslicer)
+lhinter = rve.ImagingLikelihood(interobs, skyslicer)
+# End interferometer likelihood
diff --git a/demo/mosaicing.py b/demo/mosaicing.py
index 1ddceabf..21eb5b72 100644
--- a/demo/mosaicing.py
+++ b/demo/mosaicing.py
@@ -2,13 +2,11 @@
# Copyright(C) 2019-2021 Max-Planck-Society
# Author: Philipp Arras
-import numpy as np
-from functools import reduce
-from operator import add
import matplotlib.pyplot as plt
import nifty7 as ift
import resolve as rve
+from cfg_mosaicing import lhinter, lhsd, logsky, sky, xfov, yfov
@rve.onlymaster
@@ -18,123 +16,12 @@ def imsave(fname, field):
def main():
- rve.set_nthreads(2)
- rve.set_wgridding(False)
- diffusefluxlevel = 20
-
- interobs_extended = {
- f"extended{ii}": rve.Observation.load(f"/data/mosaicing/extended_field{ii}.npz")
- for ii in range(5, 135)
- }
- interobs_compact = {
- f"compact{ii}": rve.Observation.load(f"/data/mosaicing/compact_field{ii}.npz")
- for ii in range(5, 153)
- }
- interobs = {**interobs_compact, **interobs_extended}
- sdobs = {"sd0": rve.SingleDishObservation.load("/data/mosaicing/totalpower.npz")}
- assert len(set(sdobs.keys()) & set(interobs.keys())) == 0
-
- rve.set_epsilon(1 / 10 / max(o.max_snr() for o in interobs.values()))
-
- # Compute phase center and define domain
- # Convention: The phase center is the coordinate of the [nx/2, ny/2] pixel
- xs, ys = [], []
- for vv in interobs.values():
- xs.append(vv.direction.phase_center[0])
- ys.append(vv.direction.phase_center[1])
- for vv in sdobs.values():
- foo, bar = vv.pointings.phase_centers.T
- xs.extend(foo)
- ys.extend(bar)
- xs, ys = np.array(xs), np.array(ys)
- global_phase_center_casa = -2.7180339078174175, -0.5210930907116116
- global_phase_center = (max(xs) + min(xs)) / 2, (max(ys) + min(ys)) / 2
- margin = 1 * rve.ARCMIN2RAD
- xfov = 2 * max(abs(xs - global_phase_center[0])) + 2 * margin
- yfov = 2 * max(abs(ys - global_phase_center[1])) + 2 * margin
if rve.mpi.master:
print(f"Field of view: {xfov/rve.ARCMIN2RAD:.2f}' x {yfov/rve.ARCMIN2RAD:.2f}'")
- npix = np.array([1800, 1800])
- fov = np.array([xfov, yfov])
- dom = ift.RGSpace(npix, fov / npix)
- # End compute phase center and define domain
-
- # Plot pointings
- xs, ys = [], []
- for vv in interobs.values():
- xs.append(vv.direction.phase_center[0])
- ys.append(vv.direction.phase_center[1])
- xs, ys = np.array(xs), np.array(ys)
- plt.scatter(xs, ys, label="Interferometer")
- for kk, vv in sdobs.items():
- plt.scatter(
- *vv.pointings.phase_centers.T, label="Single dish " + kk, alpha=0.2, s=1
- )
- plt.scatter(*global_phase_center_casa, label="Phase center CASA")
- plt.scatter(*global_phase_center, label="Phase center resolve")
- plt.xlim([global_phase_center[0] - xfov / 2, global_phase_center[0] + xfov / 2])
- plt.ylim([global_phase_center[1] - yfov / 2, global_phase_center[1] + yfov / 2])
- plt.gca().invert_xaxis()
- plt.gca().set_aspect("equal")
- plt.legend()
- if rve.mpi.master:
- plt.savefig("pointings.png")
- plt.close()
- # End plot pointings
-
- # Single dish likelihood
- single_dish_response = []
- for kk, oo in sdobs.items():
- assert np.min(oo.freq) / np.max(oo.freq) > 0.99
- R = rve.SingleDishResponse(
- oo,
- dom,
- lambda x: rve.alma_beam_func(10.7, 0.75, np.mean(oo.freq), x),
- global_phase_center,
- )
- single_dish_response.append(R.ducktape_left(kk))
- single_dish_response = reduce(add, single_dish_response)
- dsd = ift.MultiField.from_dict({kk: o.vis for kk, o in sdobs.items()})
- invcovsd = ift.makeOp(
- ift.MultiField.from_dict({kk: o.weight for kk, o in sdobs.items()})
- )
- lhsd = ift.GaussianEnergy(dsd, invcovsd) @ single_dish_response
- imsave("sd_dirty.png", single_dish_response.adjoint(invcovsd(dsd)))
- # End single dish likelihood
-
- logsky = ift.SimpleCorrelatedField(
- dom, diffusefluxlevel, (1, 0.1), (1.5, 1), (1.2, 0.4), (0.2, 0.2), (-1.5, 0.5)
- )
- sky = logsky.exp()
-
- # Interferometer likelihood
- beam_directions = {}
- for kk, obs in interobs.items():
- assert np.min(obs.freq) / np.max(obs.freq) > 0.99
- beam_directions[kk] = rve.BeamDirection(
- obs.direction.phase_center[0] - global_phase_center[0],
- obs.direction.phase_center[1] - global_phase_center[1],
- lambda x: rve.alma_beam_func(10.7, 0.75, np.mean(obs.freq), x),
- 0.1,
- )
- skyslicer = rve.SkySlicer(logsky.target, beam_directions)
- ift.extra.check_linear_operator(skyslicer)
- lhinter = rve.ImagingLikelihood(interobs, skyslicer)
- # End interferometer likelihood
-
lh = (lhinter + lhsd) @ sky
-
- # Plots
- R = rve.StokesIResponse(interobs, skyslicer.target) @ skyslicer
- vis = ift.MultiField.from_dict({kk: o.vis for kk, o in interobs.items()})
- weight = ift.MultiField.from_dict({kk: o.weight for kk, o in interobs.items()})
- imsave("total_dirty.png", R.adjoint(vis * weight))
- imsave("skyslicer_adjoint.png", skyslicer.adjoint(ift.full(skyslicer.target, 1.0)))
- # End plots
-
plotter = rve.Plotter("png", "plots")
plotter.add("logsky", logsky)
- plotter.add("sky", logsky.exp(), cmap="inferno_r")
+ plotter.add("sky", sky, cmap="inferno_r")
plotter.add("power spectrum logsky", logsky.power_spectrum)
# plotter.add_histogram(
# "normalized residuals (original weights)", lh.normalized_residual
diff --git a/demo/post_mosaicing.py b/demo/post_mosaicing.py
new file mode 100644
index 00000000..ca3d257a
--- /dev/null
+++ b/demo/post_mosaicing.py
@@ -0,0 +1,112 @@
+from functools import reduce
+from operator import add
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import nifty7 as ift
+import resolve as rve
+from cfg_mosaicing import (
+ global_phase_center,
+ global_phase_center_casa,
+ interobs,
+ interobs_compact,
+ interobs_extended,
+ responsesd,
+ sky,
+ sdobs,
+ skyslicer,
+ skyslicer_compact,
+ skyslicer_extended,
+ xfov,
+ yfov,
+)
+
+
+def plot_pointings():
+ xs, ys = [], []
+ for vv in interobs.values():
+ xs.append(vv.direction.phase_center[0])
+ ys.append(vv.direction.phase_center[1])
+ xs, ys = np.array(xs), np.array(ys)
+ plt.scatter(xs, ys, label="Interferometer")
+ for kk, vv in sdobs.items():
+ plt.scatter(
+ *vv.pointings.phase_centers.T, label="Single dish " + kk, alpha=0.2, s=1
+ )
+ plt.scatter(*global_phase_center_casa, label="Phase center CASA")
+ plt.scatter(*global_phase_center, label="Phase center resolve")
+ plt.xlim([global_phase_center[0] - xfov / 2, global_phase_center[0] + xfov / 2])
+ plt.ylim([global_phase_center[1] - yfov / 2, global_phase_center[1] + yfov / 2])
+ plt.gca().invert_xaxis()
+ plt.gca().set_aspect("equal")
+ plt.legend()
+ plt.savefig("pointings.png")
+ plt.close()
+
+
+def imshow(ax, fld, **kwargs):
+ kwargs = {**dict(origin="lower", cmap="inferno_r"), **kwargs}
+ return ax.imshow(fld.val.T, **kwargs)
+
+
+def imshow_symmetric(ax, fld):
+ lim = np.max(np.abs(fld.val))
+ im = imshow(ax, fld, vmin=-lim, vmax=lim, cmap="seismic")
+ plt.colorbar(im, ax=ax)
+
+
+def dirty_images(mean_sky):
+ responses = [responsesd]
+ responses.extend(
+ [
+ rve.StokesIResponse(obs, skyslc.target) @ skyslc
+ for obs, skyslc in (
+ (interobs_extended, skyslicer_extended),
+ (interobs_compact, skyslicer_compact),
+ (interobs, skyslicer),
+ )
+ ]
+ )
+ responses.append(reduce(add, responses))
+
+ data = [
+ ift.MultiField.from_dict({kk: o.vis for kk, o in obs.items()})
+ for obs in [sdobs, interobs_extended, interobs_compact, interobs]
+ ]
+ data.append(reduce(lambda a, b: a.unite(b), data))
+
+ weights = [
+ ift.MultiField.from_dict({kk: o.weight for kk, o in obs.items()})
+ for obs in [sdobs, interobs_extended, interobs_compact, interobs]
+ ]
+ weights.append(reduce(lambda a, b: a.unite(b), weights))
+
+ fig, axs = plt.subplots(5, 5, figsize=(20, 20), sharey=True, sharex=True)
+
+ for ii, (R, icov, d) in enumerate(zip(responses, weights, data)):
+ jj = 0
+ imshow_symmetric(axs[ii, jj], R.adjoint(ift.full(R.target, 1.0)))
+ jj += 1
+ imshow_symmetric(axs[ii, jj], R.adjoint(icov * d))
+ jj += 1
+ imshow_symmetric(axs[ii, jj], R.adjoint(icov.sqrt() * d))
+ jj += 1
+ residual = d - R(mean_sky)
+ imshow_symmetric(axs[ii, jj], R.adjoint(icov * residual))
+ jj += 1
+ imshow_symmetric(axs[ii, jj], R.adjoint(icov.sqrt() * residual))
+ plt.tight_layout()
+ plt.savefig("dirty_overview.png")
+
+
+def main():
+ state = rve.MinimizationState.load("iter8")
+ skysc = ift.StatCalculator()
+ for ss in state:
+ skysc.add(sky.force(ss))
+ dirty_images(skysc.mean)
+
+
+if __name__ == "__main__":
+ main()
--
GitLab
From ded45f07dee955669f5cf11153b706a6ca289790 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Wed, 10 Feb 2021 10:58:33 +0100
Subject: [PATCH 21/37] Work
---
demo/cfg_mosaicing.py | 13 +++++++++++++
demo/mosaicing.py | 31 ++++++++++++++++++++-----------
2 files changed, 33 insertions(+), 11 deletions(-)
diff --git a/demo/cfg_mosaicing.py b/demo/cfg_mosaicing.py
index d34490d4..5ac84652 100644
--- a/demo/cfg_mosaicing.py
+++ b/demo/cfg_mosaicing.py
@@ -65,6 +65,19 @@ invcovsd = ift.makeOp(
ift.MultiField.from_dict({kk: o.weight for kk, o in sdobs.items()})
)
lhsd = ift.GaussianEnergy(dsd, invcovsd) @ responsesd
+
+
+def VariableCovarianceGaussianEnergy(data, signal_response, invcov):
+ resi = ift.Adder(dsd, True).ducktape_left("resi") @ responsesd
+ dtype = {kk: vv.dtype for kk, vv in data.items()}
+ return ift.VariableCovarianceGaussianEnergy(data.domain, "resi", "icov", dtype) @ (
+ resi + invcov.ducktape_left("icov")
+ )
+
+
+assert list(dsd.domain.keys()) == ["sd0"]
+invcovop = ift.FieldAdapter(dsd.domain["sd0"], "sd0invcov").ducktape_left("sd0").exp()
+lhsd_varcov = VariableCovarianceGaussianEnergy(dsd, responsesd, invcovop)
# End single dish likelihood
# Sky model
diff --git a/demo/mosaicing.py b/demo/mosaicing.py
index 21eb5b72..cdecee56 100644
--- a/demo/mosaicing.py
+++ b/demo/mosaicing.py
@@ -4,9 +4,9 @@
import matplotlib.pyplot as plt
+import cfg_mosaicing as c
import nifty7 as ift
import resolve as rve
-from cfg_mosaicing import lhinter, lhsd, logsky, sky, xfov, yfov
@rve.onlymaster
@@ -17,15 +17,20 @@ def imsave(fname, field):
def main():
if rve.mpi.master:
- print(f"Field of view: {xfov/rve.ARCMIN2RAD:.2f}' x {yfov/rve.ARCMIN2RAD:.2f}'")
- lh = (lhinter + lhsd) @ sky
- plotter = rve.Plotter("png", "plots")
- plotter.add("logsky", logsky)
- plotter.add("sky", sky, cmap="inferno_r")
- plotter.add("power spectrum logsky", logsky.power_spectrum)
- # plotter.add_histogram(
- # "normalized residuals (original weights)", lh.normalized_residual
- # )
+ print(
+ f"Field of view: {c.xfov/rve.ARCMIN2RAD:.2f}' x {c.yfov/rve.ARCMIN2RAD:.2f}'"
+ )
+ minimize(c.lhsd @ c.sky, "sd")
+ minimize(c.lh_compact @ c.sky, "compact")
+ minimize(c.lh_extended @ c.sky, "extended")
+ minimize((c.lh_compact + c.lh_extended + c.lhsd) @ c.sky, "all")
+
+
+def minimize(lh, prefix):
+ plotter = rve.Plotter("png", f"plots_{prefix}")
+ plotter.add("logsky", c.logsky)
+ plotter.add("sky", c.sky, cmap="inferno_r")
+ plotter.add("power spectrum logsky", c.logsky.power_spectrum)
ham = ift.StandardHamiltonian(
lh, ift.AbsDeltaEnergyController(0.5, 5, 300, name="Sampling")
)
@@ -37,7 +42,11 @@ def main():
for ii in range(20):
state = rve.simple_minimize(ham, state.mean, 3, mini)
plotter.plot(f"iter{ii}", state)
- state.save(f"iter{ii}")
+ state.save(f"{prefix}iter")
+
+
+def minimize_with_varcov(lh, prefix):
+ pass
if __name__ == "__main__":
--
GitLab
From b851a9455a3c55cbaeb83ebe1220b2d1dc7e2e89 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Wed, 10 Feb 2021 13:16:57 +0100
Subject: [PATCH 22/37] Add KeyPrefixer
---
resolve/simple_operators.py | 21 +++++++++++++++++++++
test/test_general.py | 8 ++++++++
2 files changed, 29 insertions(+)
diff --git a/resolve/simple_operators.py b/resolve/simple_operators.py
index 7a900e9f..a2c294bd 100644
--- a/resolve/simple_operators.py
+++ b/resolve/simple_operators.py
@@ -64,3 +64,24 @@ class AddEmptyDimensionAtEnd(ift.LinearOperator):
else:
x = x.val[..., 0]
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()}"
diff --git a/test/test_general.py b/test/test_general.py
index 3511ce12..d09c9f2e 100644
--- a/test/test_general.py
+++ b/test/test_general.py
@@ -310,3 +310,11 @@ def test_intop():
dom = ift.RGSpace((12, 12))
op = rve.WienerIntegrations(freqdomain, dom)
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)
--
GitLab
From c3ea91888406f3db7a25a65cc485e88457eefb8e Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Wed, 10 Feb 2021 13:24:55 +0100
Subject: [PATCH 23/37] Reorder
---
demo/cfg_mosaicing.py | 14 +++++++-------
demo/mosaicing.py | 2 --
2 files changed, 7 insertions(+), 9 deletions(-)
diff --git a/demo/cfg_mosaicing.py b/demo/cfg_mosaicing.py
index 5ac84652..5a90c856 100644
--- a/demo/cfg_mosaicing.py
+++ b/demo/cfg_mosaicing.py
@@ -47,6 +47,13 @@ fov = np.array([xfov, yfov])
dom = ift.RGSpace(npix, fov / npix)
# End compute phase center and define domain
+# Sky model
+logsky = ift.SimpleCorrelatedField(
+ dom, diffusefluxlevel, (1, 0.1), (1.5, 1), (1.2, 0.4), (0.2, 0.2), (-1.5, 0.5)
+)
+sky = logsky.exp()
+# End sky model
+
# Single dish likelihood
responsesd = []
for kk, oo in sdobs.items():
@@ -80,13 +87,6 @@ invcovop = ift.FieldAdapter(dsd.domain["sd0"], "sd0invcov").ducktape_left("sd0")
lhsd_varcov = VariableCovarianceGaussianEnergy(dsd, responsesd, invcovop)
# End single dish likelihood
-# Sky model
-logsky = ift.SimpleCorrelatedField(
- dom, diffusefluxlevel, (1, 0.1), (1.5, 1), (1.2, 0.4), (0.2, 0.2), (-1.5, 0.5)
-)
-sky = logsky.exp()
-# End sky model
-
# Interferometer likelihood
def get_beamdirections(obs):
diff --git a/demo/mosaicing.py b/demo/mosaicing.py
index cdecee56..d398c66f 100644
--- a/demo/mosaicing.py
+++ b/demo/mosaicing.py
@@ -37,8 +37,6 @@ def minimize(lh, prefix):
fld = 0.1 * ift.from_random(lh.domain)
state = rve.MinimizationState(fld, [])
mini = ift.NewtonCG(ift.GradientNormController(name="newton", iteration_limit=5))
- # TODO Add minisanity to resolve
- # TODO Fix histogram plots
for ii in range(20):
state = rve.simple_minimize(ham, state.mean, 3, mini)
plotter.plot(f"iter{ii}", state)
--
GitLab
From f6e37b526fa395baa714c7537ea7f3bda2871ef5 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Wed, 10 Feb 2021 13:55:42 +0100
Subject: [PATCH 24/37] Implement variable covariance Gaussian energy
---
demo/cfg_mosaicing.py | 41 +++++++++++++++++++------------------
demo/mosaicing.py | 3 +++
demo/post_mosaicing.py | 4 ++++
demo/pre_mosaicing.py | 27 ++++++++++++++++++++++++
resolve/simple_operators.py | 24 ++++++++++++++++++++--
test/test_general.py | 3 ++-
6 files changed, 79 insertions(+), 23 deletions(-)
create mode 100644 demo/pre_mosaicing.py
diff --git a/demo/cfg_mosaicing.py b/demo/cfg_mosaicing.py
index 5a90c856..84acc7c7 100644
--- a/demo/cfg_mosaicing.py
+++ b/demo/cfg_mosaicing.py
@@ -72,21 +72,24 @@ invcovsd = ift.makeOp(
ift.MultiField.from_dict({kk: o.weight for kk, o in sdobs.items()})
)
lhsd = ift.GaussianEnergy(dsd, invcovsd) @ responsesd
-
-
-def VariableCovarianceGaussianEnergy(data, signal_response, invcov):
- resi = ift.Adder(dsd, True).ducktape_left("resi") @ responsesd
- dtype = {kk: vv.dtype for kk, vv in data.items()}
- return ift.VariableCovarianceGaussianEnergy(data.domain, "resi", "icov", dtype) @ (
- resi + invcov.ducktape_left("icov")
- )
-
-
-assert list(dsd.domain.keys()) == ["sd0"]
-invcovop = ift.FieldAdapter(dsd.domain["sd0"], "sd0invcov").ducktape_left("sd0").exp()
-lhsd_varcov = VariableCovarianceGaussianEnergy(dsd, responsesd, invcovop)
+invcovsdop = invcovsd @ rve.KeyPrefixer(dsd.domain, "invcov_inp").adjoint.exp()
+lhsd_varcov = rve.MultiDomainVariableCovarianceGaussianEnergy(
+ dsd, responsesd.ducktape("sky"), invcovsdop
+)
# End single dish likelihood
+# Consistency checks
+ift.extra.assert_allclose(
+ invcovsdop(ift.full(invcovsdop.domain, 0.0)),
+ invcovsd(ift.full(invcovsd.domain, 1.0)),
+)
+_, lhsd1 = lhsd_varcov.simplify_for_constant_input(ift.full(invcovsdop.domain, 0.0))
+foo = ift.from_random(lhsd1.domain)
+# Inaccuracy probably due to tr-log term in VariableCovGaussianEnergy
+ift.extra.assert_allclose(lhsd(foo["sky"]), lhsd1(foo), rtol=5e-3)
+del lhsd1, foo
+# End consistency check
+
# Interferometer likelihood
def get_beamdirections(obs):
@@ -102,11 +105,9 @@ def get_beamdirections(obs):
return beam_directions
-beamsc = get_beamdirections(interobs_compact)
-beamsext = get_beamdirections(interobs_extended)
-skyslicer_compact = rve.SkySlicer(logsky.target, beamsc)
-skyslicer_extended = rve.SkySlicer(logsky.target, beamsext)
-skyslicer = rve.SkySlicer(logsky.target, {**beamsc, **beamsext})
-ift.extra.check_linear_operator(skyslicer)
-lhinter = rve.ImagingLikelihood(interobs, skyslicer)
+skyslicer_compact = rve.SkySlicer(logsky.target, get_beamdirections(interobs_compact))
+skyslicer_extended = rve.SkySlicer(logsky.target, get_beamdirections(interobs_extended))
+ift.extra.check_linear_operator(skyslicer_compact)
+lh_compact = rve.ImagingLikelihood(interobs_compact, skyslicer_compact)
+lh_extended = rve.ImagingLikelihood(interobs_extended, skyslicer_extended)
# End interferometer likelihood
diff --git a/demo/mosaicing.py b/demo/mosaicing.py
index d398c66f..4776f948 100644
--- a/demo/mosaicing.py
+++ b/demo/mosaicing.py
@@ -20,6 +20,9 @@ def main():
print(
f"Field of view: {c.xfov/rve.ARCMIN2RAD:.2f}' x {c.yfov/rve.ARCMIN2RAD:.2f}'"
)
+
+ minimize_with_varcov(c.lhsd_varcov @ c.sky.ducktape_left("sky"), [], "sd")
+
minimize(c.lhsd @ c.sky, "sd")
minimize(c.lh_compact @ c.sky, "compact")
minimize(c.lh_extended @ c.sky, "extended")
diff --git a/demo/post_mosaicing.py b/demo/post_mosaicing.py
index ca3d257a..046486fc 100644
--- a/demo/post_mosaicing.py
+++ b/demo/post_mosaicing.py
@@ -1,3 +1,7 @@
+# SPDX-License-Identifier: GPL-3.0-or-later
+# Copyright(C) 2021 Max-Planck-Society
+# Author: Philipp Arras
+
from functools import reduce
from operator import add
diff --git a/demo/pre_mosaicing.py b/demo/pre_mosaicing.py
new file mode 100644
index 00000000..d9ca1d0b
--- /dev/null
+++ b/demo/pre_mosaicing.py
@@ -0,0 +1,27 @@
+# SPDX-License-Identifier: GPL-3.0-or-later
+# Copyright(C) 2021 Max-Planck-Society
+# Author: Philipp Arras
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import cfg_mosaicing as c
+import nifty7 as ift
+import resolve as rve
+
+def plot_sd_data():
+ for kk, vv in c.dsd.items():
+ lim = np.max(np.abs(vv.val))
+ plt.hist(vv.val.ravel(), bins=np.arange(-lim, lim+1, 1))
+ plt.xlabel("Flux [Jy]")
+ plt.ylabel("Count")
+ plt.title("Single dish data")
+ plt.savefig(f"single_dish_data_{kk}.png")
+ plt.close()
+
+
+def main():
+ plot_sd_data()
+
+if __name__=="__main__":
+ main()
diff --git a/resolve/simple_operators.py b/resolve/simple_operators.py
index a2c294bd..8e3b7b76 100644
--- a/resolve/simple_operators.py
+++ b/resolve/simple_operators.py
@@ -1,12 +1,15 @@
# SPDX-License-Identifier: GPL-3.0-or-later
-# Copyright(C) 2019-2020 Max-Planck-Society
+# Copyright(C) 2019-2021 Max-Planck-Society
# Author: Philipp Arras
+from functools import reduce
+from operator import add
+
import numpy as np
import nifty7 as ift
-from .util import my_assert_isinstance, my_asserteq
+from .util import my_assert, my_assert_isinstance, my_asserteq
class AddEmptyDimension(ift.LinearOperator):
@@ -85,3 +88,20 @@ class KeyPrefixer(ift.LinearOperator):
def __repr__(self):
return f"{self.domain.keys()} -> {self.target.keys()}"
+
+
+def MultiDomainVariableCovarianceGaussianEnergy(data, signal_response, invcov):
+ my_asserteq(data.domain, signal_response.target, invcov.target)
+ my_assert_isinstance(data.domain, ift.MultiDomain)
+ my_assert_isinstance(signal_response.domain, ift.MultiDomain)
+ my_assert(ift.is_operator(invcov))
+ my_assert(ift.is_operator(signal_response))
+ resi = KeyPrefixer(data.domain, "resi") @ ift.Adder(data, True) @ signal_response
+ invcov = KeyPrefixer(data.domain, "icov") @ invcov
+ res = [
+ ift.VariableCovarianceGaussianEnergy(
+ data.domain[kk], "resi" + kk, "icov" + kk, data[kk].dtype
+ )
+ for kk in data.keys()
+ ]
+ return reduce(add, res) @ (resi + invcov)
diff --git a/test/test_general.py b/test/test_general.py
index d09c9f2e..61a17224 100644
--- a/test/test_general.py
+++ b/test/test_general.py
@@ -2,8 +2,9 @@
# Copyright(C) 2020 Max-Planck-Society
# Author: Philipp Arras
-import numpy as np
from os.path import join
+
+import numpy as np
import pytest
import nifty7 as ift
--
GitLab
From d8ea77305096ea46e5088ab6520de5e90127b5be Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Wed, 10 Feb 2021 18:57:07 +0100
Subject: [PATCH 25/37] Manage implementing a somewhat working version for the
single dish data set
---
demo/cfg_mosaicing.py | 30 +++++++++++++----
demo/mosaicing.py | 28 ++++++++++++++--
demo/post_mosaicing.py | 26 +--------------
demo/pre_mosaicing.py | 65 +++++++++++++++++++++++++++++++++----
misc/convert_alma_data.py | 10 ++++++
resolve/likelihood.py | 26 +++++++++++++--
resolve/simple_operators.py | 20 ++++++++----
7 files changed, 156 insertions(+), 49 deletions(-)
diff --git a/demo/cfg_mosaicing.py b/demo/cfg_mosaicing.py
index 84acc7c7..aea45566 100644
--- a/demo/cfg_mosaicing.py
+++ b/demo/cfg_mosaicing.py
@@ -10,18 +10,20 @@ import nifty7 as ift
import resolve as rve
-diffusefluxlevel = 20
+diffusefluxlevel = 29
npix = np.array([1800, 1800])
rve.set_nthreads(2)
rve.set_wgridding(False)
interobs_extended = {
f"extended{ii}": rve.Observation.load(f"/data/mosaicing/extended_field{ii}.npz")
- for ii in range(5, 135)
+ # for ii in range(5, 135)
+ for ii in range(5, 7)
}
interobs_compact = {
f"compact{ii}": rve.Observation.load(f"/data/mosaicing/compact_field{ii}.npz")
- for ii in range(5, 153)
+ # for ii in range(5, 153)
+ for ii in range(5, 7)
}
interobs = {**interobs_compact, **interobs_extended}
sdobs = {"sd0": rve.SingleDishObservation.load("/data/mosaicing/totalpower.npz")}
@@ -49,7 +51,7 @@ dom = ift.RGSpace(npix, fov / npix)
# Sky model
logsky = ift.SimpleCorrelatedField(
- dom, diffusefluxlevel, (1, 0.1), (1.5, 1), (1.2, 0.4), (0.2, 0.2), (-1.5, 0.5)
+ dom, diffusefluxlevel, (1, 1), (1.5, 1), (1.2, 0.4), (0.2, 0.2), (-1.5, 0.5)
)
sky = logsky.exp()
# End sky model
@@ -72,7 +74,21 @@ invcovsd = ift.makeOp(
ift.MultiField.from_dict({kk: o.weight for kk, o in sdobs.items()})
)
lhsd = ift.GaussianEnergy(dsd, invcovsd) @ responsesd
-invcovsdop = invcovsd @ rve.KeyPrefixer(dsd.domain, "invcov_inp").adjoint.exp()
+
+if True:
+ # Learn one number per observation
+ foo = []
+ for kk in dsd.keys():
+ foo.append(
+ ift.ContractionOperator(invcovsd.domain[kk], (0, 1, 2))
+ .adjoint.ducktape(kk + "invcovinp")
+ .ducktape_left(kk)
+ )
+ invcovsdop = invcovsd @ reduce(add, foo).exp()
+else:
+ # Learn one number per data point
+ invcovsdop = invcovsd @ rve.KeyPrefixer(dsd.domain, "invcov_inp").adjoint.exp()
+
lhsd_varcov = rve.MultiDomainVariableCovarianceGaussianEnergy(
dsd, responsesd.ducktape("sky"), invcovsdop
)
@@ -85,8 +101,8 @@ ift.extra.assert_allclose(
)
_, lhsd1 = lhsd_varcov.simplify_for_constant_input(ift.full(invcovsdop.domain, 0.0))
foo = ift.from_random(lhsd1.domain)
-# Inaccuracy probably due to tr-log term in VariableCovGaussianEnergy
-ift.extra.assert_allclose(lhsd(foo["sky"]), lhsd1(foo), rtol=5e-3)
+# FIXME Inaccuracy probably due to tr-log term in VariableCovGaussianEnergy
+ift.extra.assert_allclose(lhsd(foo["sky"]), lhsd1(foo), rtol=1e-2)
del lhsd1, foo
# End consistency check
diff --git a/demo/mosaicing.py b/demo/mosaicing.py
index 4776f948..0285bbbf 100644
--- a/demo/mosaicing.py
+++ b/demo/mosaicing.py
@@ -21,7 +21,13 @@ def main():
f"Field of view: {c.xfov/rve.ARCMIN2RAD:.2f}' x {c.yfov/rve.ARCMIN2RAD:.2f}'"
)
- minimize_with_varcov(c.lhsd_varcov @ c.sky.ducktape_left("sky"), [], "sd")
+ # minimize(c.lhsd @ c.sky, "sd")
+ minimize_with_varcov(
+ c.lhsd_varcov.partial_insert(c.sky.ducktape_left("sky")),
+ c.invcovsdop.domain.keys(),
+ "sd",
+ )
+ exit()
minimize(c.lhsd @ c.sky, "sd")
minimize(c.lh_compact @ c.sky, "compact")
@@ -46,8 +52,24 @@ def minimize(lh, prefix):
state.save(f"{prefix}iter")
-def minimize_with_varcov(lh, prefix):
- pass
+def minimize_with_varcov(lh, varcov_keys, prefix):
+ assert set(varcov_keys) < set(lh.domain.keys())
+ plotter = rve.Plotter("png", f"plots_{prefix}")
+ plotter.add("logsky", c.logsky)
+ plotter.add("sky", c.sky, cmap="inferno_r")
+ plotter.add("power spectrum logsky", c.logsky.power_spectrum)
+ ham = ift.StandardHamiltonian(
+ lh, ift.AbsDeltaEnergyController(0.5, 5, 300, name="Sampling")
+ )
+ fld = 0.1 * ift.from_random(lh.domain)
+ state = rve.MinimizationState(fld, [])
+ mini = ift.NewtonCG(ift.GradientNormController(name="newton", iteration_limit=5))
+ for ii in range(20):
+ state = rve.simple_minimize(
+ ham, state.mean, 3, mini, constants=varcov_keys if ii < 3 else []
+ )
+ plotter.plot(f"{ii}both", state)
+ state.save(f"{prefix}both")
if __name__ == "__main__":
diff --git a/demo/post_mosaicing.py b/demo/post_mosaicing.py
index 046486fc..250ca80e 100644
--- a/demo/post_mosaicing.py
+++ b/demo/post_mosaicing.py
@@ -17,9 +17,8 @@ from cfg_mosaicing import (
interobs_compact,
interobs_extended,
responsesd,
- sky,
sdobs,
- skyslicer,
+ sky,
skyslicer_compact,
skyslicer_extended,
xfov,
@@ -27,28 +26,6 @@ from cfg_mosaicing import (
)
-def plot_pointings():
- xs, ys = [], []
- for vv in interobs.values():
- xs.append(vv.direction.phase_center[0])
- ys.append(vv.direction.phase_center[1])
- xs, ys = np.array(xs), np.array(ys)
- plt.scatter(xs, ys, label="Interferometer")
- for kk, vv in sdobs.items():
- plt.scatter(
- *vv.pointings.phase_centers.T, label="Single dish " + kk, alpha=0.2, s=1
- )
- plt.scatter(*global_phase_center_casa, label="Phase center CASA")
- plt.scatter(*global_phase_center, label="Phase center resolve")
- plt.xlim([global_phase_center[0] - xfov / 2, global_phase_center[0] + xfov / 2])
- plt.ylim([global_phase_center[1] - yfov / 2, global_phase_center[1] + yfov / 2])
- plt.gca().invert_xaxis()
- plt.gca().set_aspect("equal")
- plt.legend()
- plt.savefig("pointings.png")
- plt.close()
-
-
def imshow(ax, fld, **kwargs):
kwargs = {**dict(origin="lower", cmap="inferno_r"), **kwargs}
return ax.imshow(fld.val.T, **kwargs)
@@ -68,7 +45,6 @@ def dirty_images(mean_sky):
for obs, skyslc in (
(interobs_extended, skyslicer_extended),
(interobs_compact, skyslicer_compact),
- (interobs, skyslicer),
)
]
)
diff --git a/demo/pre_mosaicing.py b/demo/pre_mosaicing.py
index d9ca1d0b..efae1cb3 100644
--- a/demo/pre_mosaicing.py
+++ b/demo/pre_mosaicing.py
@@ -4,24 +4,77 @@
import matplotlib.pyplot as plt
import numpy as np
+import nifty7 as ift
import cfg_mosaicing as c
-import nifty7 as ift
-import resolve as rve
+
def plot_sd_data():
+ # FIXME QUESTION: In which unit is the single dish data? Jy/beam? Maybe learn relative factor between single dish and interferometer or figure out the correct factor.
for kk, vv in c.dsd.items():
- lim = np.max(np.abs(vv.val))
- plt.hist(vv.val.ravel(), bins=np.arange(-lim, lim+1, 1))
- plt.xlabel("Flux [Jy]")
+ obs = c.sdobs[kk]
+ directions = obs.pointings.phase_centers
+ ind = np.any(obs.vis.val <= 0.0, axis=(0, 2))
+ plt.scatter(*directions[ind].T, label="Non-positive", alpha=0.4, s=2)
+ plt.scatter(*directions[~ind].T, label="Positive", alpha=0.4, s=2)
+ plt.gca().invert_xaxis()
+ plt.gca().set_aspect("equal")
+ plt.legend()
+ plt.savefig(f"single_dish_negative_{kk}.png")
+ plt.close()
+
+ vv = obs.apply_flags(vv.val)
+ lim = np.max(np.abs(vv))
+ plt.hist(vv.ravel(), bins=np.arange(-lim, lim + 1, 1))
+ plt.xlabel("Flux [???]")
plt.ylabel("Count")
plt.title("Single dish data")
plt.savefig(f"single_dish_data_{kk}.png")
plt.close()
+def plot_pointings():
+ xs, ys = [], []
+ for vv in c.interobs.values():
+ xs.append(vv.direction.phase_center[0])
+ ys.append(vv.direction.phase_center[1])
+ xs, ys = np.array(xs), np.array(ys)
+ plt.scatter(xs, ys, label="Interferometer")
+ for kk, vv in c.sdobs.items():
+ plt.scatter(
+ *vv.pointings.phase_centers.T, label="Single dish " + kk, alpha=0.2, s=1
+ )
+ plt.scatter(*c.global_phase_center_casa, label="Phase center CASA")
+ plt.scatter(*c.global_phase_center, label="Phase center resolve")
+ plt.xlim(
+ [c.global_phase_center[0] - c.xfov / 2, c.global_phase_center[0] + c.xfov / 2]
+ )
+ plt.ylim(
+ [c.global_phase_center[1] - c.yfov / 2, c.global_phase_center[1] + c.yfov / 2]
+ )
+ plt.gca().invert_xaxis()
+ plt.gca().set_aspect("equal")
+ plt.legend()
+ plt.savefig("pointings.png")
+ plt.close()
+
+
+def plot_sky_prior():
+ p = ift.Plot()
+ pspecs = []
+ for _ in range(8):
+ pos = ift.from_random(c.sky.domain)
+ p.add(c.logsky(pos))
+ pspecs.append(c.logsky.power_spectrum.force(pos))
+ p.add(pspecs)
+ p.output(name="sky_prior.png")
+
+
def main():
+ plot_sky_prior()
+ plot_pointings()
plot_sd_data()
-if __name__=="__main__":
+
+if __name__ == "__main__":
main()
diff --git a/misc/convert_alma_data.py b/misc/convert_alma_data.py
index 0962f982..4a638632 100644
--- a/misc/convert_alma_data.py
+++ b/misc/convert_alma_data.py
@@ -28,6 +28,7 @@ def cleanup(weight, vis, freq, uvw=None):
mywgt = mywgt[:, ind]
myvis = myvis[:, ind]
myfreq = freq[ind]
+ myvis[mywgt == 0.0] = 0.0
if uvw is None:
return mywgt, myvis, myfreq
return mywgt, myvis, myfreq, myuvw
@@ -81,6 +82,15 @@ for f, spectral_window, single_dish, fieldid in [
del flag
vis = sel("FLOAT_DATA" if single_dish else "DATA")[:, channel_slice]
+ # FIXME QUESTION Why is there negative data?
+ if single_dish:
+ assert not np.iscomplexobj(vis)
+ ind = vis * np.sqrt(wgt) <= 2.0
+ print(
+ f"{np.sum(ind)} data points are not 2-sigma positive ({np.sum(ind)/ind.size*100:.1f} %)"
+ )
+ wgt[ind] = 0.0
+
# Average polarization
assert vis.shape[2] == 2
assert wgt.shape[2] == 2
diff --git a/resolve/likelihood.py b/resolve/likelihood.py
index c714f35c..a659bbbe 100644
--- a/resolve/likelihood.py
+++ b/resolve/likelihood.py
@@ -2,6 +2,9 @@
# Copyright(C) 2020 Max-Planck-Society
# Author: Philipp Arras
+from functools import reduce
+from operator import add
+
import numpy as np
import nifty7 as ift
@@ -14,7 +17,6 @@ from .util import my_assert, my_assert_isinstance, my_asserteq
def _get_mask(observation):
# Only needed for variable covariance gaussian energy
my_assert_isinstance(observation, Observation)
-
vis = observation.vis
flags = observation.flags
if not np.any(flags):
@@ -23,6 +25,24 @@ def _get_mask(observation):
return mask, mask(vis), mask(observation.weight)
+def get_mask_multi_field(weight):
+ assert isinstance(weight, ift.MultiField)
+ op = []
+ for kk, ww in weight.items():
+ flags = ww.val == 0.0
+ if np.any(flags):
+ op.append(
+ ift.MaskOperator(ift.makeField(ww.domain, flags))
+ .ducktape(kk)
+ .ducktape_left(kk)
+ )
+ else:
+ op.append(ift.ScalingOperator(ww.domain, 1.0))
+ op = reduce(add, op)
+ assert op.domain == weight.domain
+ return op
+
+
def _Likelihood(operator, normalized_residual_operator):
my_assert_isinstance(operator, ift.Operator)
my_asserteq(operator.target, ift.DomainTuple.scalar_domain())
@@ -160,7 +180,9 @@ def ImagingLikelihood(
if inverse_covariance_operator is None:
if mosaicing:
vis = ift.MultiField.from_dict({kk: o.vis for kk, o in observation.items()})
- weight = ift.MultiField.from_dict({kk: o.weight for kk, o in observation.items()})
+ weight = ift.MultiField.from_dict(
+ {kk: o.weight for kk, o in observation.items()}
+ )
else:
vis, weight = observation.vis, observation.weight
return _build_gauss_lh_nres(model_data, vis, weight)
diff --git a/resolve/simple_operators.py b/resolve/simple_operators.py
index 8e3b7b76..ad834f4c 100644
--- a/resolve/simple_operators.py
+++ b/resolve/simple_operators.py
@@ -91,17 +91,25 @@ class KeyPrefixer(ift.LinearOperator):
def MultiDomainVariableCovarianceGaussianEnergy(data, signal_response, invcov):
+ from .likelihood import get_mask_multi_field
+
my_asserteq(data.domain, signal_response.target, invcov.target)
my_assert_isinstance(data.domain, ift.MultiDomain)
my_assert_isinstance(signal_response.domain, ift.MultiDomain)
my_assert(ift.is_operator(invcov))
my_assert(ift.is_operator(signal_response))
+ res = []
+ invcovfld = invcov(ift.full(invcov.domain, 1.0))
+ mask = get_mask_multi_field(invcovfld)
+ data = mask(data)
+ signal_response = mask @ signal_response
+ invcov = mask @ invcov
+ for kk in data.keys():
+ res.append(
+ ift.VariableCovarianceGaussianEnergy(
+ 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
- res = [
- ift.VariableCovarianceGaussianEnergy(
- data.domain[kk], "resi" + kk, "icov" + kk, data[kk].dtype
- )
- for kk in data.keys()
- ]
return reduce(add, res) @ (resi + invcov)
--
GitLab
From 60f4c40a84d837b6bf25550925bbbc564393e85a Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Fri, 12 Feb 2021 11:07:12 +0100
Subject: [PATCH 26/37] Add debugging script
---
demo/debug_mosaicing.py | 51 +++++++++++++++++++++++++++++++++++++++++
1 file changed, 51 insertions(+)
create mode 100644 demo/debug_mosaicing.py
diff --git a/demo/debug_mosaicing.py b/demo/debug_mosaicing.py
new file mode 100644
index 00000000..eec4303a
--- /dev/null
+++ b/demo/debug_mosaicing.py
@@ -0,0 +1,51 @@
+import resolve as rve
+import numpy as np
+import matplotlib.pyplot as plt
+import nifty7 as ift
+import cfg_mosaicing as c
+
+
+if False:
+ state = rve.MinimizationState.load("sditer")
+ bins = np.linspace(-20, 20, 40)
+ for kk in c.invcovsd.target.keys():
+ for ss in state:
+ arr = (
+ (c.invcovsd.get_sqrt() @ ift.Adder(c.dsd, True) @ c.responsesd @ c.sky)
+ .force(ss)
+ .val
+ )
+ arr = c.sdobs[kk].apply_flags(arr[kk]).ravel()
+ plt.hist(arr.clip(-20, 20), bins=bins, alpha=0.5)
+ print("Reduced chi**2:", np.sum(arr ** 2) / arr.size)
+ plt.title(f"Noise weighted residual {kk}")
+ plt.savefig(f"debug_nwresi_{kk}.png")
+ plt.close()
+else:
+ state = rve.MinimizationState.load("sdboth")
+
+ upperlim = 10
+ bins = np.linspace(0, upperlim, upperlim * 10)
+ for kk in c.invcovsdop.target.keys():
+ for ss in state:
+ learnedsd = 1 / c.invcovsdop.force(ss)[kk].sqrt()
+ datasd = 1 / c.invcovsd(ift.full(c.invcovsd.domain, 1.0))[kk].sqrt()
+ plt.hist(
+ (learnedsd / datasd).sqrt().val.ravel().clip(0, upperlim), bins=bins
+ )
+ plt.axvline(1.0)
+ plt.xlabel("Learned sigma / sigma in measurement set")
+ plt.savefig(f"debug{kk}.png")
+ plt.close()
+
+ bins = np.linspace(-20, 20, 40)
+ for kk in c.invcovsd.target.keys():
+ for ss in state:
+ arr = (ift.Adder(c.dsd, True) @ c.responsesd @ c.sky).force(ss)
+ arr = arr * c.invcovsdop.force(ss).sqrt()
+ arr = c.sdobs[kk].apply_flags(arr.val[kk]).ravel()
+ plt.hist(arr.clip(-20, 20), bins=bins, alpha=0.5)
+ print("Reduced chi**2:", np.sum(arr ** 2) / arr.size)
+ plt.title(f"Noise weighted residual {kk}")
+ plt.savefig(f"debug_nwresi_{kk}.png")
+ plt.close()
--
GitLab
From 5a85042f90f6722557031e36e30b0ee8d5910ae5 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Fri, 12 Feb 2021 13:47:22 +0100
Subject: [PATCH 27/37] Cleanup
---
resolve/likelihood.py | 30 ------------------------------
1 file changed, 30 deletions(-)
diff --git a/resolve/likelihood.py b/resolve/likelihood.py
index a659bbbe..859e4477 100644
--- a/resolve/likelihood.py
+++ b/resolve/likelihood.py
@@ -83,36 +83,6 @@ def _varcov(observation, Rs, inverse_covariance_operator):
return _build_varcov_gauss_lh_nres(residual, inverse_covariance_operator, dtype)
-def SingleDishLikelihood(
- observation,
- sky_operator,
- inverse_covariance_operator=None,
- calibration_operator=None,
-):
- if inverse_covariance_operator is not None or calibration_operator is not None:
- raise NotImplementedError
- mosaicing = isinstance(observation, dict)
- if mosaicing:
- vis = ift.MultiField.from_dict({kk: o.vis for kk, o in observation.items()})
- weight = ift.MultiField.from_dict(
- {kk: o.weight for kk, o in observation.items()}
- )
- single_dish_response = []
- for kk, oo in observation.items():
- assert np.min(oo.freq) / np.max(oo.freq) > 0.99
- R = rve.SingleDishResponse(
- oo,
- sky_operator.target,
- lambda x: rve.alma_beam_func(10.7, 0.75, np.mean(oo.freq), x),
- global_phase_center,
- )
- single_dish_response.append(R.ducktape_left(kk))
- from functools import reduce
- from operator import add
-
- single_dish_response = reduce(add, single_dish_response)
-
-
def ImagingLikelihood(
observation,
sky_operator,
--
GitLab
From 14e299d70b70070b1d6d677f98696334e2463a87 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Fri, 12 Feb 2021 15:17:01 +0100
Subject: [PATCH 28/37] Sky target always MultiDomain
---
demo/cfg_mosaicing.py | 12 +++++++-----
resolve/likelihood.py | 9 +++------
resolve/mosaicing/single_dish_response.py | 2 +-
3 files changed, 11 insertions(+), 12 deletions(-)
diff --git a/demo/cfg_mosaicing.py b/demo/cfg_mosaicing.py
index aea45566..fc1713e9 100644
--- a/demo/cfg_mosaicing.py
+++ b/demo/cfg_mosaicing.py
@@ -53,7 +53,7 @@ dom = ift.RGSpace(npix, fov / npix)
logsky = ift.SimpleCorrelatedField(
dom, diffusefluxlevel, (1, 1), (1.5, 1), (1.2, 0.4), (0.2, 0.2), (-1.5, 0.5)
)
-sky = logsky.exp()
+sky = logsky.exp().ducktape_left("sky")
# End sky model
# Single dish likelihood
@@ -90,7 +90,7 @@ else:
invcovsdop = invcovsd @ rve.KeyPrefixer(dsd.domain, "invcov_inp").adjoint.exp()
lhsd_varcov = rve.MultiDomainVariableCovarianceGaussianEnergy(
- dsd, responsesd.ducktape("sky"), invcovsdop
+ dsd, responsesd, invcovsdop
)
# End single dish likelihood
@@ -102,7 +102,7 @@ ift.extra.assert_allclose(
_, lhsd1 = lhsd_varcov.simplify_for_constant_input(ift.full(invcovsdop.domain, 0.0))
foo = ift.from_random(lhsd1.domain)
# FIXME Inaccuracy probably due to tr-log term in VariableCovGaussianEnergy
-ift.extra.assert_allclose(lhsd(foo["sky"]), lhsd1(foo), rtol=1e-2)
+ift.extra.assert_allclose(lhsd(foo), lhsd1(foo), rtol=1e-2)
del lhsd1, foo
# End consistency check
@@ -124,6 +124,8 @@ def get_beamdirections(obs):
skyslicer_compact = rve.SkySlicer(logsky.target, get_beamdirections(interobs_compact))
skyslicer_extended = rve.SkySlicer(logsky.target, get_beamdirections(interobs_extended))
ift.extra.check_linear_operator(skyslicer_compact)
-lh_compact = rve.ImagingLikelihood(interobs_compact, skyslicer_compact)
-lh_extended = rve.ImagingLikelihood(interobs_extended, skyslicer_extended)
+lh_compact = rve.ImagingLikelihood(interobs_compact, skyslicer_compact).ducktape("sky")
+lh_extended = rve.ImagingLikelihood(interobs_extended, skyslicer_extended).ducktape(
+ "sky"
+)
# End interferometer likelihood
diff --git a/resolve/likelihood.py b/resolve/likelihood.py
index 859e4477..49d1b4b1 100644
--- a/resolve/likelihood.py
+++ b/resolve/likelihood.py
@@ -31,13 +31,10 @@ def get_mask_multi_field(weight):
for kk, ww in weight.items():
flags = ww.val == 0.0
if np.any(flags):
- op.append(
- ift.MaskOperator(ift.makeField(ww.domain, flags))
- .ducktape(kk)
- .ducktape_left(kk)
- )
+ myop = ift.MaskOperator(ift.makeField(ww.domain, flags))
else:
- op.append(ift.ScalingOperator(ww.domain, 1.0))
+ myop = ift.ScalingOperator(ww.domain, 1.0)
+ op.append(myop.ducktape(kk).ducktape_left(kk))
op = reduce(add, op)
assert op.domain == weight.domain
return op
diff --git a/resolve/mosaicing/single_dish_response.py b/resolve/mosaicing/single_dish_response.py
index d4aa2b87..e54f0660 100644
--- a/resolve/mosaicing/single_dish_response.py
+++ b/resolve/mosaicing/single_dish_response.py
@@ -21,7 +21,7 @@ def SingleDishResponse(observation, domain, beam_function, global_phase_center):
fld = ift.from_random(conv.domain)
ift.extra.assert_allclose(conv(fld).integrate(), fld.integrate())
pc = observation.pointings.phase_centers.T - np.array(global_phase_center)[:, None]
- pc = pc + (np.array(domain.shape)*np.array(domain[0].distances)/2)[:, None]
+ pc = pc + (np.array(domain.shape) * np.array(domain[0].distances) / 2)[:, None]
# Convention: pointing convention (see also BeamDirection)
pc[0] *= -1
interp = ift.LinearInterpolator(domain, pc)
--
GitLab
From 374d3d36d1bcca5211645dbbe128338b009ac92e Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Fri, 12 Feb 2021 15:17:32 +0100
Subject: [PATCH 29/37] Add additive term for dealing with negative data
---
demo/cfg_mosaicing.py | 4 ++++
demo/mosaicing.py | 3 +++
resolve/mosaicing/single_dish_response.py | 11 +++++++++--
3 files changed, 16 insertions(+), 2 deletions(-)
diff --git a/demo/cfg_mosaicing.py b/demo/cfg_mosaicing.py
index fc1713e9..a9ae0af9 100644
--- a/demo/cfg_mosaicing.py
+++ b/demo/cfg_mosaicing.py
@@ -58,6 +58,9 @@ sky = logsky.exp().ducktape_left("sky")
# Single dish likelihood
responsesd = []
+additive_term = ift.SimpleCorrelatedField(
+ dom, 0, (1, 1), (1, 1), (1, 1), None, (-2, 0.5), "negativity_madness"
+)
for kk, oo in sdobs.items():
# Single dish response does not support multifrequency
assert np.min(oo.freq) / np.max(oo.freq) > 0.99
@@ -66,6 +69,7 @@ for kk, oo in sdobs.items():
dom,
lambda x: rve.alma_beam_func(10.7, 0.75, np.mean(oo.freq), x),
global_phase_center,
+ additive_term=additive_term,
)
responsesd.append(R.ducktape_left(kk))
responsesd = reduce(add, responsesd)
diff --git a/demo/mosaicing.py b/demo/mosaicing.py
index 0285bbbf..8e6484ac 100644
--- a/demo/mosaicing.py
+++ b/demo/mosaicing.py
@@ -40,6 +40,9 @@ def minimize(lh, prefix):
plotter.add("logsky", c.logsky)
plotter.add("sky", c.sky, cmap="inferno_r")
plotter.add("power spectrum logsky", c.logsky.power_spectrum)
+ plotter.add("additive term", c.additive_term)
+ plotter.add("additive term power spectrum", c.additive_term.power_spectrum)
+
ham = ift.StandardHamiltonian(
lh, ift.AbsDeltaEnergyController(0.5, 5, 300, name="Sampling")
)
diff --git a/resolve/mosaicing/single_dish_response.py b/resolve/mosaicing/single_dish_response.py
index e54f0660..1d584882 100644
--- a/resolve/mosaicing/single_dish_response.py
+++ b/resolve/mosaicing/single_dish_response.py
@@ -9,7 +9,9 @@ import nifty7 as ift
from ..observation import SingleDishObservation
-def SingleDishResponse(observation, domain, beam_function, global_phase_center):
+def SingleDishResponse(
+ observation, domain, beam_function, global_phase_center, additive_term=None
+):
assert isinstance(observation, SingleDishObservation)
domain = ift.makeDomain(domain)
assert len(domain) == 1
@@ -20,6 +22,7 @@ def SingleDishResponse(observation, domain, beam_function, global_phase_center):
# FIXME Move into tests
fld = ift.from_random(conv.domain)
ift.extra.assert_allclose(conv(fld).integrate(), fld.integrate())
+
pc = observation.pointings.phase_centers.T - np.array(global_phase_center)[:, None]
pc = pc + (np.array(domain.shape) * np.array(domain[0].distances) / 2)[:, None]
# Convention: pointing convention (see also BeamDirection)
@@ -29,4 +32,8 @@ def SingleDishResponse(observation, domain, beam_function, global_phase_center):
# NOTE The volume factor above `domain.total_volume()` and the volume factor
# below `domain[0].scalar_dvol` cancel each other. They are left in the
# code such that the convolution leaves the integral invariant.
- return bc @ interp @ conv.scale(domain[0].scalar_dvol)
+
+ convsky = conv.scale(domain[0].scalar_dvol).ducktape("sky")
+ if additive_term is not None:
+ convsky = convsky + additive_term
+ return bc @ interp @ convsky
--
GitLab
From bebdc5cb2d7bd2dcd37bdd610810ca2f2c286856 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Fri, 12 Feb 2021 17:24:23 +0100
Subject: [PATCH 30/37] Use all single dish data again
---
misc/convert_alma_data.py | 16 ++++++++--------
1 file changed, 8 insertions(+), 8 deletions(-)
diff --git a/misc/convert_alma_data.py b/misc/convert_alma_data.py
index 4a638632..fd578ae2 100644
--- a/misc/convert_alma_data.py
+++ b/misc/convert_alma_data.py
@@ -82,14 +82,14 @@ for f, spectral_window, single_dish, fieldid in [
del flag
vis = sel("FLOAT_DATA" if single_dish else "DATA")[:, channel_slice]
- # FIXME QUESTION Why is there negative data?
- if single_dish:
- assert not np.iscomplexobj(vis)
- ind = vis * np.sqrt(wgt) <= 2.0
- print(
- f"{np.sum(ind)} data points are not 2-sigma positive ({np.sum(ind)/ind.size*100:.1f} %)"
- )
- wgt[ind] = 0.0
+ # # FIXME QUESTION Why is there negative data?
+ # if single_dish:
+ # assert not np.iscomplexobj(vis)
+ # ind = vis * np.sqrt(wgt) <= 2.0
+ # print(
+ # f"{np.sum(ind)} data points are not 2-sigma positive ({np.sum(ind)/ind.size*100:.1f} %)"
+ # )
+ # wgt[ind] = 0.0
# Average polarization
assert vis.shape[2] == 2
--
GitLab
From 03fe871b84195aeebc39f16f64274b996b94ebc9 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Fri, 12 Feb 2021 17:24:31 +0100
Subject: [PATCH 31/37] Kind of working version for single dish
---
demo/debug_mosaicing.py | 45 +++++++++++++++++++++++++++--------------
demo/mosaicing.py | 17 +++++++++-------
2 files changed, 40 insertions(+), 22 deletions(-)
diff --git a/demo/debug_mosaicing.py b/demo/debug_mosaicing.py
index eec4303a..1c561bc8 100644
--- a/demo/debug_mosaicing.py
+++ b/demo/debug_mosaicing.py
@@ -22,27 +22,42 @@ if False:
plt.savefig(f"debug_nwresi_{kk}.png")
plt.close()
else:
+ state = rve.MinimizationState.load("sditer")
state = rve.MinimizationState.load("sdboth")
- upperlim = 10
- bins = np.linspace(0, upperlim, upperlim * 10)
- for kk in c.invcovsdop.target.keys():
- for ss in state:
- learnedsd = 1 / c.invcovsdop.force(ss)[kk].sqrt()
- datasd = 1 / c.invcovsd(ift.full(c.invcovsd.domain, 1.0))[kk].sqrt()
- plt.hist(
- (learnedsd / datasd).sqrt().val.ravel().clip(0, upperlim), bins=bins
- )
- plt.axvline(1.0)
- plt.xlabel("Learned sigma / sigma in measurement set")
- plt.savefig(f"debug{kk}.png")
- plt.close()
+ plotter = rve.Plotter("png", "debug")
+ plotter.add("additive term", c.additive_term)
+ plotter.add("sky", ift.FieldAdapter(c.sky.target["sky"], "sky")@ c.sky, cmap="inferno_r")
+ plotter.add("additive term power spectrum", c.additive_term.power_spectrum)
+ plotter.plot("debug", state)
+
+ learnnoise = True
+ if learnnoise:
+ upperlim = 10
+ bins = np.linspace(0, upperlim, upperlim * 10)
+ for kk in c.invcovsdop.target.keys():
+ for ss in state:
+ learnedsd = 1 / c.invcovsdop.force(ss)[kk].sqrt()
+ datasd = 1 / c.invcovsd(ift.full(c.invcovsd.domain, 1.0))[kk].sqrt()
+ plt.hist(
+ (learnedsd / datasd).sqrt().val.ravel().clip(0, upperlim), bins=bins
+ )
+ plt.axvline(1.0)
+ plt.xlabel("Learned sigma / sigma in measurement set")
+ plt.savefig(f"debug{kk}.png")
+ plt.close()
+
bins = np.linspace(-20, 20, 40)
for kk in c.invcovsd.target.keys():
for ss in state:
- arr = (ift.Adder(c.dsd, True) @ c.responsesd @ c.sky).force(ss)
- arr = arr * c.invcovsdop.force(ss).sqrt()
+ arr = (ift.Adder(c.dsd, True) @ c.responsesd.partial_insert(c.sky)).force(
+ ss
+ )
+ if learnnoise:
+ arr = arr * c.invcovsdop.force(ss).sqrt()
+ else:
+ arr = c.invcovsd.get_sqrt()(arr)
arr = c.sdobs[kk].apply_flags(arr.val[kk]).ravel()
plt.hist(arr.clip(-20, 20), bins=bins, alpha=0.5)
print("Reduced chi**2:", np.sum(arr ** 2) / arr.size)
diff --git a/demo/mosaicing.py b/demo/mosaicing.py
index 8e6484ac..61151979 100644
--- a/demo/mosaicing.py
+++ b/demo/mosaicing.py
@@ -21,15 +21,12 @@ def main():
f"Field of view: {c.xfov/rve.ARCMIN2RAD:.2f}' x {c.yfov/rve.ARCMIN2RAD:.2f}'"
)
- # minimize(c.lhsd @ c.sky, "sd")
+ # minimize(c.lhsd.partial_insert(c.sky), "sd")
minimize_with_varcov(
- c.lhsd_varcov.partial_insert(c.sky.ducktape_left("sky")),
- c.invcovsdop.domain.keys(),
- "sd",
+ c.lhsd_varcov.partial_insert(c.sky), c.invcovsdop.domain.keys(), "sd"
)
exit()
- minimize(c.lhsd @ c.sky, "sd")
minimize(c.lh_compact @ c.sky, "compact")
minimize(c.lh_extended @ c.sky, "extended")
minimize((c.lh_compact + c.lh_extended + c.lhsd) @ c.sky, "all")
@@ -38,7 +35,9 @@ def main():
def minimize(lh, prefix):
plotter = rve.Plotter("png", f"plots_{prefix}")
plotter.add("logsky", c.logsky)
- plotter.add("sky", c.sky, cmap="inferno_r")
+ plotter.add(
+ "sky", ift.FieldAdapter(c.sky.target["sky"], "sky") @ c.sky, cmap="inferno_r"
+ )
plotter.add("power spectrum logsky", c.logsky.power_spectrum)
plotter.add("additive term", c.additive_term)
plotter.add("additive term power spectrum", c.additive_term.power_spectrum)
@@ -59,8 +58,12 @@ def minimize_with_varcov(lh, varcov_keys, prefix):
assert set(varcov_keys) < set(lh.domain.keys())
plotter = rve.Plotter("png", f"plots_{prefix}")
plotter.add("logsky", c.logsky)
- plotter.add("sky", c.sky, cmap="inferno_r")
+ plotter.add(
+ "sky", ift.FieldAdapter(c.sky.target["sky"], "sky") @ c.sky, cmap="inferno_r"
+ )
plotter.add("power spectrum logsky", c.logsky.power_spectrum)
+ plotter.add("additive term", c.additive_term)
+ plotter.add("additive term power spectrum", c.additive_term.power_spectrum)
ham = ift.StandardHamiltonian(
lh, ift.AbsDeltaEnergyController(0.5, 5, 300, name="Sampling")
)
--
GitLab
From 2d65c643af85ea4bb0a7dc3d5bf47a4f86777287 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Mon, 15 Feb 2021 14:35:08 +0100
Subject: [PATCH 32/37] Add comparison to ALMA reconstruction
---
demo/cfg_mosaicing.py | 3 +-
demo/post_mosaicing.py | 64 +++++++++++++++++---------
resolve/__init__.py | 2 +-
resolve/fits.py | 100 ++++++++++++++++++++++++++++-------------
4 files changed, 113 insertions(+), 56 deletions(-)
diff --git a/demo/cfg_mosaicing.py b/demo/cfg_mosaicing.py
index a9ae0af9..ec512ef8 100644
--- a/demo/cfg_mosaicing.py
+++ b/demo/cfg_mosaicing.py
@@ -53,7 +53,8 @@ dom = ift.RGSpace(npix, fov / npix)
logsky = ift.SimpleCorrelatedField(
dom, diffusefluxlevel, (1, 1), (1.5, 1), (1.2, 0.4), (0.2, 0.2), (-1.5, 0.5)
)
-sky = logsky.exp().ducktape_left("sky")
+sky_domaintuple = logsky.exp()
+sky = sky_domaintuple.ducktape_left("sky")
# End sky model
# Single dish likelihood
diff --git a/demo/post_mosaicing.py b/demo/post_mosaicing.py
index 250ca80e..67a2c01a 100644
--- a/demo/post_mosaicing.py
+++ b/demo/post_mosaicing.py
@@ -10,20 +10,8 @@ import numpy as np
import nifty7 as ift
import resolve as rve
-from cfg_mosaicing import (
- global_phase_center,
- global_phase_center_casa,
- interobs,
- interobs_compact,
- interobs_extended,
- responsesd,
- sdobs,
- sky,
- skyslicer_compact,
- skyslicer_extended,
- xfov,
- yfov,
-)
+
+import cfg_mosaicing as c
def imshow(ax, fld, **kwargs):
@@ -31,6 +19,12 @@ def imshow(ax, fld, **kwargs):
return ax.imshow(fld.val.T, **kwargs)
+def imshow_with_colorbar(ax, fld, **kwargs):
+ kwargs = {**dict(origin="lower", cmap="inferno_r"), **kwargs}
+ im = ax.imshow(fld.val.T, **kwargs)
+ plt.colorbar(im, ax=ax)
+
+
def imshow_symmetric(ax, fld):
lim = np.max(np.abs(fld.val))
im = imshow(ax, fld, vmin=-lim, vmax=lim, cmap="seismic")
@@ -38,13 +32,13 @@ def imshow_symmetric(ax, fld):
def dirty_images(mean_sky):
- responses = [responsesd]
+ responses = [c.responsesd]
responses.extend(
[
- rve.StokesIResponse(obs, skyslc.target) @ skyslc
+ rve.StokesIResponse(obs, skyslc.target) @ skyslc.ducktape("sky")
for obs, skyslc in (
- (interobs_extended, skyslicer_extended),
- (interobs_compact, skyslicer_compact),
+ (c.interobs_extended, c.skyslicer_extended),
+ (c.interobs_compact, c.skyslicer_compact),
)
]
)
@@ -52,13 +46,13 @@ def dirty_images(mean_sky):
data = [
ift.MultiField.from_dict({kk: o.vis for kk, o in obs.items()})
- for obs in [sdobs, interobs_extended, interobs_compact, interobs]
+ for obs in [c.sdobs, c.interobs_extended, c.interobs_compact, c.interobs]
]
data.append(reduce(lambda a, b: a.unite(b), data))
weights = [
ift.MultiField.from_dict({kk: o.weight for kk, o in obs.items()})
- for obs in [sdobs, interobs_extended, interobs_compact, interobs]
+ for obs in [c.sdobs, c.interobs_extended, c.interobs_compact, c.interobs]
]
weights.append(reduce(lambda a, b: a.unite(b), weights))
@@ -81,12 +75,38 @@ def dirty_images(mean_sky):
def main():
- state = rve.MinimizationState.load("iter8")
+ state = rve.MinimizationState.load("sdboth")
skysc = ift.StatCalculator()
for ss in state:
- skysc.add(sky.force(ss))
+ skysc.add(c.sky_domaintuple.force(ss))
+
+ sd_comparison(skysc.mean)
dirty_images(skysc.mean)
+def sd_comparison(mean_sky):
+ fld, refval = rve.fits2field("../resources_multipointing/M83.spw17.I.sd.fits")
+ # FIXME Check with convert_alma_data.py script
+ minfreq = 230172900000.0
+ ind = int(np.round((minfreq - refval[0]) / fld.domain[0].distances[0]))
+ extractor = ift.DomainTupleFieldInserter(fld.domain, 0, (ind,)).adjoint
+ casasky = extractor(fld)
+
+ fig, axs = plt.subplots(1, 2)
+ axs = list(axs)
+
+ axx = axs.pop(0)
+ imshow_with_colorbar(axx, casasky)
+ axx.set_title("CASA single dish")
+
+ axx = axs.pop(0)
+ imshow_with_colorbar(axx, mean_sky)
+ axx.set_title("Resolve single dish posterior mean")
+
+ plt.tight_layout()
+ fig.savefig("sd_reconstructions.png")
+ plt.close()
+
+
if __name__ == "__main__":
main()
diff --git a/resolve/__init__.py b/resolve/__init__.py
index b74b7fda..d69bc4b0 100644
--- a/resolve/__init__.py
+++ b/resolve/__init__.py
@@ -2,7 +2,7 @@ from .antenna_positions import AntennaPositions
from .calibration import calibration_distribution
from .constants import *
from .direction import *
-from .fits import field2fits
+from .fits import field2fits, fits2field
from .global_config import *
from .likelihood import *
from .minimization import Minimization, MinimizationState, simple_minimize
diff --git a/resolve/fits.py b/resolve/fits.py
index 4f2ecc6d..93700159 100644
--- a/resolve/fits.py
+++ b/resolve/fits.py
@@ -1,5 +1,5 @@
# SPDX-License-Identifier: GPL-3.0-or-later
-# Copyright(C) 2019-2020 Max-Planck-Society
+# Copyright(C) 2019-2021 Max-Planck-Society
# Author: Philipp Arras
import time
@@ -7,6 +7,9 @@ from os.path import splitext
import numpy as np
+import nifty7 as ift
+
+from .constants import DEG2RAD
from .mpi import onlymaster
@@ -41,35 +44,68 @@ def field2fits(field, file_name, overwrite, direction=None):
base, ext = splitext(file_name)
hdulist.writeto(base + ext, overwrite=overwrite)
- # @staticmethod
- # def make_from_file(file_name):
- # with pyfits.open(file_name) as hdu_list:
- # lst = hdu_list[0]
- # pcx = lst.header['CRVAL1']/180*np.pi
- # pcy = lst.header['CRVAL2']/180*np.pi
- # equ = lst.header['EQUINOX']
- # return FitsWriter([pcx, pcy], equ)
- # @staticmethod
- # def fits2field(file_name, ignore_units=False, from_wsclean=False):
- # with pyfits.open(file_name) as hdu_list:
- # image_data = np.squeeze(hdu_list[0].data).astype(np.float64)
- # head = hdu_list[0].header
- # dstx = abs(head['CDELT1']*np.pi/180)
- # dsty = abs(head['CDELT2']*np.pi/180)
- # if not ignore_units:
- # if head['BUNIT'] == 'JY/BEAM':
- # fac = np.pi/4/np.log(2)
- # scale = fac*head['BMAJ']*head['BMIN']*(np.pi/180)**2
- # elif head['BUNIT'] == 'JY/PIXEL':
- # scale = dstx*dsty
- # else:
- # scale = 1
- # image_data /= scale
- # if from_wsclean:
- # image_data = image_data[::-1].T[:, :-1]
- # image_data = np.pad(image_data, ((0, 0), (1, 0)), mode='constant')
- # else:
- # image_data = image_data.T[:, ::-1]
- # dom = ift.RGSpace(image_data.shape, (dstx, dsty))
- # return ift.makeField(dom, image_data)
+def fits2field(file_name, ignore_units=False, from_wsclean=False):
+ import astropy.io.fits as pyfits
+
+ with pyfits.open(file_name) as hdu_list:
+ image_data = hdu_list[0].data.astype(np.float64)
+ assert image_data.shape[0] == 1 # Only one Stokes component
+ image_data = image_data[0]
+ head = hdu_list[0].header
+ assert head["CUNIT1"].strip() == "deg"
+ assert head["CUNIT2"].strip() == "deg"
+ assert head["CUNIT3"].strip() == "Hz"
+ refs = []
+ refs.append([float(head["CRVAL3"]), int(head["CRPIX3"]), head["CDELT3"]])
+ refs.append(
+ [
+ float(head["CRVAL2"]) * DEG2RAD,
+ int(head["CRPIX2"]),
+ head["CDELT2"] * DEG2RAD,
+ ]
+ )
+ refs.append(
+ [
+ float(head["CRVAL1"]) * DEG2RAD,
+ int(head["CRPIX1"]),
+ head["CDELT1"] * DEG2RAD,
+ ]
+ )
+
+ if not ignore_units:
+ if head["BUNIT"].upper() == "JY/BEAM":
+ fac = np.pi / 4 / np.log(2)
+ scale = fac * head["BMAJ"] * head["BMIN"] * (np.pi / 180) ** 2
+ elif head["BUNIT"].upper() == "JY/PIXEL":
+ scale = abs(refs[0][2] * refs[1][2])
+ else:
+ scale = 1
+ image_data /= scale
+
+ # Convert CASA conventions to resolve conventions
+ inds = 0, 2, 1
+ image_data = np.transpose(image_data, inds)
+ refs = [refs[ii] for ii in inds]
+ refs[1][2] *= -1
+
+ for ii, (_, mypx, mydst) in enumerate(refs):
+ if mydst == 0.0:
+ raise RuntimeError
+ if mydst > 0:
+ continue
+ image_data = np.flip(image_data, ii)
+ refs[ii][2] *= -1
+ # FIXME Assume pixel counting start at 0. Maybe also 1?
+ refs[ii][1] = image_data.shape[ii] - mypx
+
+ refval = tuple(refs[ii][0] for ii in range(3))
+ refpx = tuple(refs[ii][1] for ii in range(3))
+ dsts = tuple(refs[ii][2] for ii in range(3))
+ dom = (
+ ift.RGSpace(image_data.shape[0], dsts[0]),
+ ift.RGSpace(image_data.shape[1:], dsts[1:]),
+ )
+ refval = tuple(refval[ii] - refpx[ii] * dsts[ii] for ii in range(3))
+ del refpx
+ return ift.makeField(dom, image_data), refval
--
GitLab
From e4b6373f757173a211e4021fbbf61ad1171ea0af Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Mon, 15 Feb 2021 14:35:28 +0100
Subject: [PATCH 33/37] Add caching for primary beam
---
demo/cfg_mosaicing.py | 6 ++++--
resolve/primary_beam.py | 17 +++++++++++++++--
2 files changed, 19 insertions(+), 4 deletions(-)
diff --git a/demo/cfg_mosaicing.py b/demo/cfg_mosaicing.py
index ec512ef8..11b2013f 100644
--- a/demo/cfg_mosaicing.py
+++ b/demo/cfg_mosaicing.py
@@ -68,7 +68,7 @@ for kk, oo in sdobs.items():
R = rve.SingleDishResponse(
oo,
dom,
- lambda x: rve.alma_beam_func(10.7, 0.75, np.mean(oo.freq), x),
+ lambda x: rve.alma_beam_func(10.7, 0.75, np.mean(oo.freq), x, use_cache=True),
global_phase_center,
additive_term=additive_term,
)
@@ -120,7 +120,9 @@ def get_beamdirections(obs):
beam_directions[kk] = rve.BeamDirection(
oo.direction.phase_center[0] - global_phase_center[0],
oo.direction.phase_center[1] - global_phase_center[1],
- lambda x: rve.alma_beam_func(10.7, 0.75, np.mean(oo.freq), x),
+ lambda x: rve.alma_beam_func(
+ 10.7, 0.75, np.mean(oo.freq), x, use_cache=True
+ ),
0.1,
)
return beam_directions
diff --git a/resolve/primary_beam.py b/resolve/primary_beam.py
index f627f5e2..e13e7738 100644
--- a/resolve/primary_beam.py
+++ b/resolve/primary_beam.py
@@ -179,11 +179,24 @@ def vla_beam(domain, freq):
return ift.makeOp(ift.makeField(dom, beam))
-def alma_beam_func(D, d, freq, x):
+def alma_beam_func(D, d, freq, x, use_cache=False):
assert isinstance(x, np.ndarray)
assert x.ndim < 3
assert np.max(np.abs(x)) < np.pi / np.sqrt(2)
+ if not use_cache:
+ return _compute_alma_beam(D, d, freq, x)
+
+ iden = "_".join([str(ll) for ll in [D, d, freq]] + [str(ll) for ll in x.shape])
+ fname = f".beamcache{iden}.npy"
+ try:
+ return np.load(fname)
+ except FileNotFoundError:
+ arr = _compute_alma_beam(D, d, freq, x)
+ np.save(fname, arr)
+
+
+def _compute_alma_beam(D, d, freq, x):
a = freq / SPEEDOFLIGHT
b = d / D
x = np.pi * a * D * x
@@ -191,4 +204,4 @@ def alma_beam_func(D, d, freq, x):
x[mask] = 1
sol = 2 / (x * (1 - b ** 2)) * (sc.jn(1, x) - b * sc.jn(1, x * b))
sol[mask] = 1
- return sol ** 2
+ return sol * sol
--
GitLab
From 68aaabba29318d6a88eb549a40dd7f224b971873 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Mon, 15 Feb 2021 14:35:41 +0100
Subject: [PATCH 34/37] Speed up
---
demo/cfg_mosaicing.py | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/demo/cfg_mosaicing.py b/demo/cfg_mosaicing.py
index 11b2013f..78e7c073 100644
--- a/demo/cfg_mosaicing.py
+++ b/demo/cfg_mosaicing.py
@@ -18,12 +18,12 @@ rve.set_wgridding(False)
interobs_extended = {
f"extended{ii}": rve.Observation.load(f"/data/mosaicing/extended_field{ii}.npz")
# for ii in range(5, 135)
- for ii in range(5, 7)
+ for ii in range(5, 6)
}
interobs_compact = {
f"compact{ii}": rve.Observation.load(f"/data/mosaicing/compact_field{ii}.npz")
# for ii in range(5, 153)
- for ii in range(5, 7)
+ for ii in range(5, 6)
}
interobs = {**interobs_compact, **interobs_extended}
sdobs = {"sd0": rve.SingleDishObservation.load("/data/mosaicing/totalpower.npz")}
--
GitLab
From 63ce5d0c986f2e13943151cc53547e199a063e1b Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Tue, 16 Feb 2021 22:22:50 +0100
Subject: [PATCH 35/37] Work on plots
---
demo/post_mosaicing.py | 78 +++++++++++++++++++++++++++++++++++++-----
1 file changed, 70 insertions(+), 8 deletions(-)
diff --git a/demo/post_mosaicing.py b/demo/post_mosaicing.py
index 67a2c01a..ead57523 100644
--- a/demo/post_mosaicing.py
+++ b/demo/post_mosaicing.py
@@ -80,28 +80,90 @@ def main():
for ss in state:
skysc.add(c.sky_domaintuple.force(ss))
- sd_comparison(skysc.mean)
+ logskysc = ift.StatCalculator()
+ for ss in state:
+ logskysc.add(c.logsky.force(ss))
+
+ sc_sd_simdata = ift.StatCalculator()
+ op = c.responsesd.partial_insert(c.sky)
+ for ss in state:
+ sc_sd_simdata.add(op.force(ss))
+ sd_comparison(skysc.mean, logskysc.mean, sc_sd_simdata)
dirty_images(skysc.mean)
-def sd_comparison(mean_sky):
+def sd_comparison(mean_sky, mean_logsky, sc_sd_simdata):
fld, refval = rve.fits2field("../resources_multipointing/M83.spw17.I.sd.fits")
- # FIXME Check with convert_alma_data.py script
minfreq = 230172900000.0
ind = int(np.round((minfreq - refval[0]) / fld.domain[0].distances[0]))
extractor = ift.DomainTupleFieldInserter(fld.domain, 0, (ind,)).adjoint
casasky = extractor(fld)
- fig, axs = plt.subplots(1, 2)
- axs = list(axs)
+ obs = c.sdobs["sd0"]
+ r_casa = (
+ rve.SingleDishResponse(
+ obs,
+ casasky.domain,
+ lambda x: rve.alma_beam_func(10.7, 0.75, np.mean(obs.freq), x),
+ c.global_phase_center,
+ )
+ @ ift.FieldAdapter(casasky.domain, "sky").adjoint
+ )
+
+ tmp = casasky.val_rw()
+ tmp[np.isnan(tmp)] = 0.0
+ tmp = ift.makeField(casasky.domain, tmp)
+
+ # FIXME Question: How can I compute data residuals for single dish data?
+ # FIXME Question: How do I arrive at the correct units?
+ d0 = obs.vis.val.ravel()
+ d1 = r_casa(tmp).val.ravel()
+ scaling_factor = np.vdot(d0, d1) / np.vdot(d1, d1)
+ print(
+ "Crazy scaling factor to make CASA reconstruction fit the data better",
+ scaling_factor,
+ )
+
+ tmp = tmp * scaling_factor
+ d1 = r_casa(tmp).val.ravel()
+
+ plt.hist(d0)
+ plt.savefig("data.png")
+ plt.close()
+
+ plt.hist(d1)
+ plt.savefig("data_sim.png")
+ plt.close()
+
+ fig, axs = plt.subplots(2, 2)
+ axs = list(axs.ravel())
+
+ axx = axs.pop(0)
+ imshow_symmetric(axx, r_casa.adjoint(obs.vis))
+ axx.set_title("Data")
axx = axs.pop(0)
- imshow_with_colorbar(axx, casasky)
+ imshow_symmetric(axx, r_casa.adjoint(obs.vis - r_casa(tmp)))
+ axx.set_title("CASA Residual")
+
+ plt.tight_layout()
+ plt.savefig("casa_sd_reconstruction.png")
+ plt.close()
+
+ fig, axs = plt.subplots(2, 2)
+ axs = list(axs.ravel())
+
+ axx = axs.pop(0)
+ imshow_with_colorbar(axx, casasky * scaling_factor, cmap="inferno")
axx.set_title("CASA single dish")
axx = axs.pop(0)
- imshow_with_colorbar(axx, mean_sky)
- axx.set_title("Resolve single dish posterior mean")
+ imshow_with_colorbar(axx, mean_sky, cmap="inferno")
+ axx.set_title("Resolve single dish")
+
+ axx = axs.pop(0)
+ imshow_with_colorbar(axx, mean_logsky.exp(), cmap="inferno")
+ axx.set_title("Resolve log. single dish")
plt.tight_layout()
fig.savefig("sd_reconstructions.png")
--
GitLab
From 991a592aa8f0b8524695d8b574b5aa436cb77b42 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Tue, 16 Feb 2021 22:23:04 +0100
Subject: [PATCH 36/37] Various tweaks
---
demo/mosaicing.py | 3 ++-
demo/post_mosaicing.py | 13 ++++++++++++-
2 files changed, 14 insertions(+), 2 deletions(-)
diff --git a/demo/mosaicing.py b/demo/mosaicing.py
index 61151979..9ffb7387 100644
--- a/demo/mosaicing.py
+++ b/demo/mosaicing.py
@@ -68,11 +68,12 @@ def minimize_with_varcov(lh, varcov_keys, prefix):
lh, ift.AbsDeltaEnergyController(0.5, 5, 300, name="Sampling")
)
fld = 0.1 * ift.from_random(lh.domain)
+ fld = ift.MultiField.union([fld, fld.extract_by_keys(varcov_keys) * 0.0])
state = rve.MinimizationState(fld, [])
mini = ift.NewtonCG(ift.GradientNormController(name="newton", iteration_limit=5))
for ii in range(20):
state = rve.simple_minimize(
- ham, state.mean, 3, mini, constants=varcov_keys if ii < 3 else []
+ ham, state.mean, 5, mini, constants=varcov_keys if ii < 2 else []
)
plotter.plot(f"{ii}both", state)
state.save(f"{prefix}both")
diff --git a/demo/post_mosaicing.py b/demo/post_mosaicing.py
index ead57523..1567860c 100644
--- a/demo/post_mosaicing.py
+++ b/demo/post_mosaicing.py
@@ -15,17 +15,23 @@ import cfg_mosaicing as c
def imshow(ax, fld, **kwargs):
+ if isinstance(fld.domain, ift.MultiDomain):
+ fld = ift.FieldAdapter(fld.domain, "sky")(fld)
kwargs = {**dict(origin="lower", cmap="inferno_r"), **kwargs}
return ax.imshow(fld.val.T, **kwargs)
def imshow_with_colorbar(ax, fld, **kwargs):
+ if isinstance(fld.domain, ift.MultiDomain):
+ fld = ift.FieldAdapter(fld.domain, "sky")(fld)
kwargs = {**dict(origin="lower", cmap="inferno_r"), **kwargs}
- im = ax.imshow(fld.val.T, **kwargs)
+ im = ax.imshow(fld.T, **kwargs)
plt.colorbar(im, ax=ax)
def imshow_symmetric(ax, fld):
+ if isinstance(fld.domain, ift.MultiDomain):
+ fld = ift.FieldAdapter(fld.domain, "sky")(fld)
lim = np.max(np.abs(fld.val))
im = imshow(ax, fld, vmin=-lim, vmax=lim, cmap="seismic")
plt.colorbar(im, ax=ax)
@@ -56,6 +62,11 @@ def dirty_images(mean_sky):
]
weights.append(reduce(lambda a, b: a.unite(b), weights))
+ # Remove single dish response since it is not linear if additive_term is taken into account
+ responses = responses[1:]
+ data = data[1:]
+ weights = weights[1:]
+
fig, axs = plt.subplots(5, 5, figsize=(20, 20), sharey=True, sharex=True)
for ii, (R, icov, d) in enumerate(zip(responses, weights, data)):
--
GitLab
From e9d8c0a702d63fa4fcd646396eb9410b8b5c298e Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Sat, 12 Jun 2021 13:10:27 +0200
Subject: [PATCH 37/37] Cleanup
---
demo/cfg_mosaicing.py | 138 ----------------------------
demo/debug_mosaicing.py | 66 --------------
demo/mosaicing.py | 83 -----------------
demo/post_mosaicing.py | 185 --------------------------------------
demo/pre_mosaicing.py | 80 -----------------
misc/convert_alma_data.py | 141 -----------------------------
6 files changed, 693 deletions(-)
delete mode 100644 demo/cfg_mosaicing.py
delete mode 100644 demo/debug_mosaicing.py
delete mode 100644 demo/mosaicing.py
delete mode 100644 demo/post_mosaicing.py
delete mode 100644 demo/pre_mosaicing.py
delete mode 100644 misc/convert_alma_data.py
diff --git a/demo/cfg_mosaicing.py b/demo/cfg_mosaicing.py
deleted file mode 100644
index 78e7c073..00000000
--- a/demo/cfg_mosaicing.py
+++ /dev/null
@@ -1,138 +0,0 @@
-# SPDX-License-Identifier: GPL-3.0-or-later
-# Copyright(C) 2019-2021 Max-Planck-Society
-# Author: Philipp Arras
-
-import numpy as np
-from functools import reduce
-from operator import add
-
-import nifty7 as ift
-import resolve as rve
-
-
-diffusefluxlevel = 29
-npix = np.array([1800, 1800])
-
-rve.set_nthreads(2)
-rve.set_wgridding(False)
-interobs_extended = {
- f"extended{ii}": rve.Observation.load(f"/data/mosaicing/extended_field{ii}.npz")
- # for ii in range(5, 135)
- for ii in range(5, 6)
-}
-interobs_compact = {
- f"compact{ii}": rve.Observation.load(f"/data/mosaicing/compact_field{ii}.npz")
- # for ii in range(5, 153)
- for ii in range(5, 6)
-}
-interobs = {**interobs_compact, **interobs_extended}
-sdobs = {"sd0": rve.SingleDishObservation.load("/data/mosaicing/totalpower.npz")}
-assert len(set(sdobs.keys()) & set(interobs.keys())) == 0
-rve.set_epsilon(1 / 10 / max(o.max_snr() for o in interobs.values()))
-# Compute phase center and define domain
-# Convention: The phase center is the coordinate of the [nx/2, ny/2] pixel
-xs, ys = [], []
-for vv in interobs.values():
- xs.append(vv.direction.phase_center[0])
- ys.append(vv.direction.phase_center[1])
-for vv in sdobs.values():
- foo, bar = vv.pointings.phase_centers.T
- xs.extend(foo)
- ys.extend(bar)
-xs, ys = np.array(xs), np.array(ys)
-global_phase_center_casa = -2.7180339078174175, -0.5210930907116116
-global_phase_center = (max(xs) + min(xs)) / 2, (max(ys) + min(ys)) / 2
-margin = 1 * rve.ARCMIN2RAD
-xfov = 2 * max(abs(xs - global_phase_center[0])) + 2 * margin
-yfov = 2 * max(abs(ys - global_phase_center[1])) + 2 * margin
-fov = np.array([xfov, yfov])
-dom = ift.RGSpace(npix, fov / npix)
-# End compute phase center and define domain
-
-# Sky model
-logsky = ift.SimpleCorrelatedField(
- dom, diffusefluxlevel, (1, 1), (1.5, 1), (1.2, 0.4), (0.2, 0.2), (-1.5, 0.5)
-)
-sky_domaintuple = logsky.exp()
-sky = sky_domaintuple.ducktape_left("sky")
-# End sky model
-
-# Single dish likelihood
-responsesd = []
-additive_term = ift.SimpleCorrelatedField(
- dom, 0, (1, 1), (1, 1), (1, 1), None, (-2, 0.5), "negativity_madness"
-)
-for kk, oo in sdobs.items():
- # Single dish response does not support multifrequency
- assert np.min(oo.freq) / np.max(oo.freq) > 0.99
- R = rve.SingleDishResponse(
- oo,
- dom,
- lambda x: rve.alma_beam_func(10.7, 0.75, np.mean(oo.freq), x, use_cache=True),
- global_phase_center,
- additive_term=additive_term,
- )
- responsesd.append(R.ducktape_left(kk))
-responsesd = reduce(add, responsesd)
-dsd = ift.MultiField.from_dict({kk: o.vis for kk, o in sdobs.items()})
-invcovsd = ift.makeOp(
- ift.MultiField.from_dict({kk: o.weight for kk, o in sdobs.items()})
-)
-lhsd = ift.GaussianEnergy(dsd, invcovsd) @ responsesd
-
-if True:
- # Learn one number per observation
- foo = []
- for kk in dsd.keys():
- foo.append(
- ift.ContractionOperator(invcovsd.domain[kk], (0, 1, 2))
- .adjoint.ducktape(kk + "invcovinp")
- .ducktape_left(kk)
- )
- invcovsdop = invcovsd @ reduce(add, foo).exp()
-else:
- # Learn one number per data point
- invcovsdop = invcovsd @ rve.KeyPrefixer(dsd.domain, "invcov_inp").adjoint.exp()
-
-lhsd_varcov = rve.MultiDomainVariableCovarianceGaussianEnergy(
- dsd, responsesd, invcovsdop
-)
-# End single dish likelihood
-
-# Consistency checks
-ift.extra.assert_allclose(
- invcovsdop(ift.full(invcovsdop.domain, 0.0)),
- invcovsd(ift.full(invcovsd.domain, 1.0)),
-)
-_, lhsd1 = lhsd_varcov.simplify_for_constant_input(ift.full(invcovsdop.domain, 0.0))
-foo = ift.from_random(lhsd1.domain)
-# FIXME Inaccuracy probably due to tr-log term in VariableCovGaussianEnergy
-ift.extra.assert_allclose(lhsd(foo), lhsd1(foo), rtol=1e-2)
-del lhsd1, foo
-# End consistency check
-
-
-# Interferometer likelihood
-def get_beamdirections(obs):
- beam_directions = {}
- for kk, oo in obs.items():
- assert np.min(oo.freq) / np.max(oo.freq) > 0.99
- beam_directions[kk] = rve.BeamDirection(
- oo.direction.phase_center[0] - global_phase_center[0],
- oo.direction.phase_center[1] - global_phase_center[1],
- lambda x: rve.alma_beam_func(
- 10.7, 0.75, np.mean(oo.freq), x, use_cache=True
- ),
- 0.1,
- )
- return beam_directions
-
-
-skyslicer_compact = rve.SkySlicer(logsky.target, get_beamdirections(interobs_compact))
-skyslicer_extended = rve.SkySlicer(logsky.target, get_beamdirections(interobs_extended))
-ift.extra.check_linear_operator(skyslicer_compact)
-lh_compact = rve.ImagingLikelihood(interobs_compact, skyslicer_compact).ducktape("sky")
-lh_extended = rve.ImagingLikelihood(interobs_extended, skyslicer_extended).ducktape(
- "sky"
-)
-# End interferometer likelihood
diff --git a/demo/debug_mosaicing.py b/demo/debug_mosaicing.py
deleted file mode 100644
index 1c561bc8..00000000
--- a/demo/debug_mosaicing.py
+++ /dev/null
@@ -1,66 +0,0 @@
-import resolve as rve
-import numpy as np
-import matplotlib.pyplot as plt
-import nifty7 as ift
-import cfg_mosaicing as c
-
-
-if False:
- state = rve.MinimizationState.load("sditer")
- bins = np.linspace(-20, 20, 40)
- for kk in c.invcovsd.target.keys():
- for ss in state:
- arr = (
- (c.invcovsd.get_sqrt() @ ift.Adder(c.dsd, True) @ c.responsesd @ c.sky)
- .force(ss)
- .val
- )
- arr = c.sdobs[kk].apply_flags(arr[kk]).ravel()
- plt.hist(arr.clip(-20, 20), bins=bins, alpha=0.5)
- print("Reduced chi**2:", np.sum(arr ** 2) / arr.size)
- plt.title(f"Noise weighted residual {kk}")
- plt.savefig(f"debug_nwresi_{kk}.png")
- plt.close()
-else:
- state = rve.MinimizationState.load("sditer")
- state = rve.MinimizationState.load("sdboth")
-
- plotter = rve.Plotter("png", "debug")
- plotter.add("additive term", c.additive_term)
- plotter.add("sky", ift.FieldAdapter(c.sky.target["sky"], "sky")@ c.sky, cmap="inferno_r")
- plotter.add("additive term power spectrum", c.additive_term.power_spectrum)
- plotter.plot("debug", state)
-
- learnnoise = True
- if learnnoise:
- upperlim = 10
- bins = np.linspace(0, upperlim, upperlim * 10)
- for kk in c.invcovsdop.target.keys():
- for ss in state:
- learnedsd = 1 / c.invcovsdop.force(ss)[kk].sqrt()
- datasd = 1 / c.invcovsd(ift.full(c.invcovsd.domain, 1.0))[kk].sqrt()
- plt.hist(
- (learnedsd / datasd).sqrt().val.ravel().clip(0, upperlim), bins=bins
- )
- plt.axvline(1.0)
- plt.xlabel("Learned sigma / sigma in measurement set")
- plt.savefig(f"debug{kk}.png")
- plt.close()
-
-
- bins = np.linspace(-20, 20, 40)
- for kk in c.invcovsd.target.keys():
- for ss in state:
- arr = (ift.Adder(c.dsd, True) @ c.responsesd.partial_insert(c.sky)).force(
- ss
- )
- if learnnoise:
- arr = arr * c.invcovsdop.force(ss).sqrt()
- else:
- arr = c.invcovsd.get_sqrt()(arr)
- arr = c.sdobs[kk].apply_flags(arr.val[kk]).ravel()
- plt.hist(arr.clip(-20, 20), bins=bins, alpha=0.5)
- print("Reduced chi**2:", np.sum(arr ** 2) / arr.size)
- plt.title(f"Noise weighted residual {kk}")
- plt.savefig(f"debug_nwresi_{kk}.png")
- plt.close()
diff --git a/demo/mosaicing.py b/demo/mosaicing.py
deleted file mode 100644
index 9ffb7387..00000000
--- a/demo/mosaicing.py
+++ /dev/null
@@ -1,83 +0,0 @@
-# SPDX-License-Identifier: GPL-3.0-or-later
-# Copyright(C) 2019-2021 Max-Planck-Society
-# Author: Philipp Arras
-
-import matplotlib.pyplot as plt
-
-import cfg_mosaicing as c
-import nifty7 as ift
-import resolve as rve
-
-
-@rve.onlymaster
-def imsave(fname, field):
- plt.imsave(fname, field.val.T, origin="lower")
- plt.close()
-
-
-def main():
- if rve.mpi.master:
- print(
- f"Field of view: {c.xfov/rve.ARCMIN2RAD:.2f}' x {c.yfov/rve.ARCMIN2RAD:.2f}'"
- )
-
- # minimize(c.lhsd.partial_insert(c.sky), "sd")
- minimize_with_varcov(
- c.lhsd_varcov.partial_insert(c.sky), c.invcovsdop.domain.keys(), "sd"
- )
- exit()
-
- minimize(c.lh_compact @ c.sky, "compact")
- minimize(c.lh_extended @ c.sky, "extended")
- minimize((c.lh_compact + c.lh_extended + c.lhsd) @ c.sky, "all")
-
-
-def minimize(lh, prefix):
- plotter = rve.Plotter("png", f"plots_{prefix}")
- plotter.add("logsky", c.logsky)
- plotter.add(
- "sky", ift.FieldAdapter(c.sky.target["sky"], "sky") @ c.sky, cmap="inferno_r"
- )
- plotter.add("power spectrum logsky", c.logsky.power_spectrum)
- plotter.add("additive term", c.additive_term)
- plotter.add("additive term power spectrum", c.additive_term.power_spectrum)
-
- ham = ift.StandardHamiltonian(
- lh, ift.AbsDeltaEnergyController(0.5, 5, 300, name="Sampling")
- )
- fld = 0.1 * ift.from_random(lh.domain)
- state = rve.MinimizationState(fld, [])
- mini = ift.NewtonCG(ift.GradientNormController(name="newton", iteration_limit=5))
- for ii in range(20):
- state = rve.simple_minimize(ham, state.mean, 3, mini)
- plotter.plot(f"iter{ii}", state)
- state.save(f"{prefix}iter")
-
-
-def minimize_with_varcov(lh, varcov_keys, prefix):
- assert set(varcov_keys) < set(lh.domain.keys())
- plotter = rve.Plotter("png", f"plots_{prefix}")
- plotter.add("logsky", c.logsky)
- plotter.add(
- "sky", ift.FieldAdapter(c.sky.target["sky"], "sky") @ c.sky, cmap="inferno_r"
- )
- plotter.add("power spectrum logsky", c.logsky.power_spectrum)
- plotter.add("additive term", c.additive_term)
- plotter.add("additive term power spectrum", c.additive_term.power_spectrum)
- ham = ift.StandardHamiltonian(
- lh, ift.AbsDeltaEnergyController(0.5, 5, 300, name="Sampling")
- )
- fld = 0.1 * ift.from_random(lh.domain)
- fld = ift.MultiField.union([fld, fld.extract_by_keys(varcov_keys) * 0.0])
- state = rve.MinimizationState(fld, [])
- mini = ift.NewtonCG(ift.GradientNormController(name="newton", iteration_limit=5))
- for ii in range(20):
- state = rve.simple_minimize(
- ham, state.mean, 5, mini, constants=varcov_keys if ii < 2 else []
- )
- plotter.plot(f"{ii}both", state)
- state.save(f"{prefix}both")
-
-
-if __name__ == "__main__":
- main()
diff --git a/demo/post_mosaicing.py b/demo/post_mosaicing.py
deleted file mode 100644
index 1567860c..00000000
--- a/demo/post_mosaicing.py
+++ /dev/null
@@ -1,185 +0,0 @@
-# SPDX-License-Identifier: GPL-3.0-or-later
-# Copyright(C) 2021 Max-Planck-Society
-# Author: Philipp Arras
-
-from functools import reduce
-from operator import add
-
-import matplotlib.pyplot as plt
-import numpy as np
-
-import nifty7 as ift
-import resolve as rve
-
-import cfg_mosaicing as c
-
-
-def imshow(ax, fld, **kwargs):
- if isinstance(fld.domain, ift.MultiDomain):
- fld = ift.FieldAdapter(fld.domain, "sky")(fld)
- kwargs = {**dict(origin="lower", cmap="inferno_r"), **kwargs}
- return ax.imshow(fld.val.T, **kwargs)
-
-
-def imshow_with_colorbar(ax, fld, **kwargs):
- if isinstance(fld.domain, ift.MultiDomain):
- fld = ift.FieldAdapter(fld.domain, "sky")(fld)
- kwargs = {**dict(origin="lower", cmap="inferno_r"), **kwargs}
- im = ax.imshow(fld.T, **kwargs)
- plt.colorbar(im, ax=ax)
-
-
-def imshow_symmetric(ax, fld):
- if isinstance(fld.domain, ift.MultiDomain):
- fld = ift.FieldAdapter(fld.domain, "sky")(fld)
- lim = np.max(np.abs(fld.val))
- im = imshow(ax, fld, vmin=-lim, vmax=lim, cmap="seismic")
- plt.colorbar(im, ax=ax)
-
-
-def dirty_images(mean_sky):
- responses = [c.responsesd]
- responses.extend(
- [
- rve.StokesIResponse(obs, skyslc.target) @ skyslc.ducktape("sky")
- for obs, skyslc in (
- (c.interobs_extended, c.skyslicer_extended),
- (c.interobs_compact, c.skyslicer_compact),
- )
- ]
- )
- responses.append(reduce(add, responses))
-
- data = [
- ift.MultiField.from_dict({kk: o.vis for kk, o in obs.items()})
- for obs in [c.sdobs, c.interobs_extended, c.interobs_compact, c.interobs]
- ]
- data.append(reduce(lambda a, b: a.unite(b), data))
-
- weights = [
- ift.MultiField.from_dict({kk: o.weight for kk, o in obs.items()})
- for obs in [c.sdobs, c.interobs_extended, c.interobs_compact, c.interobs]
- ]
- weights.append(reduce(lambda a, b: a.unite(b), weights))
-
- # Remove single dish response since it is not linear if additive_term is taken into account
- responses = responses[1:]
- data = data[1:]
- weights = weights[1:]
-
- fig, axs = plt.subplots(5, 5, figsize=(20, 20), sharey=True, sharex=True)
-
- for ii, (R, icov, d) in enumerate(zip(responses, weights, data)):
- jj = 0
- imshow_symmetric(axs[ii, jj], R.adjoint(ift.full(R.target, 1.0)))
- jj += 1
- imshow_symmetric(axs[ii, jj], R.adjoint(icov * d))
- jj += 1
- imshow_symmetric(axs[ii, jj], R.adjoint(icov.sqrt() * d))
- jj += 1
- residual = d - R(mean_sky)
- imshow_symmetric(axs[ii, jj], R.adjoint(icov * residual))
- jj += 1
- imshow_symmetric(axs[ii, jj], R.adjoint(icov.sqrt() * residual))
- plt.tight_layout()
- plt.savefig("dirty_overview.png")
-
-
-def main():
- state = rve.MinimizationState.load("sdboth")
- skysc = ift.StatCalculator()
- for ss in state:
- skysc.add(c.sky_domaintuple.force(ss))
-
- logskysc = ift.StatCalculator()
- for ss in state:
- logskysc.add(c.logsky.force(ss))
-
- sc_sd_simdata = ift.StatCalculator()
- op = c.responsesd.partial_insert(c.sky)
- for ss in state:
- sc_sd_simdata.add(op.force(ss))
- sd_comparison(skysc.mean, logskysc.mean, sc_sd_simdata)
- dirty_images(skysc.mean)
-
-
-def sd_comparison(mean_sky, mean_logsky, sc_sd_simdata):
- fld, refval = rve.fits2field("../resources_multipointing/M83.spw17.I.sd.fits")
- minfreq = 230172900000.0
- ind = int(np.round((minfreq - refval[0]) / fld.domain[0].distances[0]))
- extractor = ift.DomainTupleFieldInserter(fld.domain, 0, (ind,)).adjoint
- casasky = extractor(fld)
-
- obs = c.sdobs["sd0"]
- r_casa = (
- rve.SingleDishResponse(
- obs,
- casasky.domain,
- lambda x: rve.alma_beam_func(10.7, 0.75, np.mean(obs.freq), x),
- c.global_phase_center,
- )
- @ ift.FieldAdapter(casasky.domain, "sky").adjoint
- )
-
- tmp = casasky.val_rw()
- tmp[np.isnan(tmp)] = 0.0
- tmp = ift.makeField(casasky.domain, tmp)
-
- # FIXME Question: How can I compute data residuals for single dish data?
- # FIXME Question: How do I arrive at the correct units?
- d0 = obs.vis.val.ravel()
- d1 = r_casa(tmp).val.ravel()
- scaling_factor = np.vdot(d0, d1) / np.vdot(d1, d1)
- print(
- "Crazy scaling factor to make CASA reconstruction fit the data better",
- scaling_factor,
- )
-
- tmp = tmp * scaling_factor
- d1 = r_casa(tmp).val.ravel()
-
- plt.hist(d0)
- plt.savefig("data.png")
- plt.close()
-
- plt.hist(d1)
- plt.savefig("data_sim.png")
- plt.close()
-
- fig, axs = plt.subplots(2, 2)
- axs = list(axs.ravel())
-
- axx = axs.pop(0)
- imshow_symmetric(axx, r_casa.adjoint(obs.vis))
- axx.set_title("Data")
-
- axx = axs.pop(0)
- imshow_symmetric(axx, r_casa.adjoint(obs.vis - r_casa(tmp)))
- axx.set_title("CASA Residual")
-
- plt.tight_layout()
- plt.savefig("casa_sd_reconstruction.png")
- plt.close()
-
- fig, axs = plt.subplots(2, 2)
- axs = list(axs.ravel())
-
- axx = axs.pop(0)
- imshow_with_colorbar(axx, casasky * scaling_factor, cmap="inferno")
- axx.set_title("CASA single dish")
-
- axx = axs.pop(0)
- imshow_with_colorbar(axx, mean_sky, cmap="inferno")
- axx.set_title("Resolve single dish")
-
- axx = axs.pop(0)
- imshow_with_colorbar(axx, mean_logsky.exp(), cmap="inferno")
- axx.set_title("Resolve log. single dish")
-
- plt.tight_layout()
- fig.savefig("sd_reconstructions.png")
- plt.close()
-
-
-if __name__ == "__main__":
- main()
diff --git a/demo/pre_mosaicing.py b/demo/pre_mosaicing.py
deleted file mode 100644
index efae1cb3..00000000
--- a/demo/pre_mosaicing.py
+++ /dev/null
@@ -1,80 +0,0 @@
-# SPDX-License-Identifier: GPL-3.0-or-later
-# Copyright(C) 2021 Max-Planck-Society
-# Author: Philipp Arras
-
-import matplotlib.pyplot as plt
-import numpy as np
-import nifty7 as ift
-
-import cfg_mosaicing as c
-
-
-def plot_sd_data():
- # FIXME QUESTION: In which unit is the single dish data? Jy/beam? Maybe learn relative factor between single dish and interferometer or figure out the correct factor.
- for kk, vv in c.dsd.items():
- obs = c.sdobs[kk]
- directions = obs.pointings.phase_centers
- ind = np.any(obs.vis.val <= 0.0, axis=(0, 2))
- plt.scatter(*directions[ind].T, label="Non-positive", alpha=0.4, s=2)
- plt.scatter(*directions[~ind].T, label="Positive", alpha=0.4, s=2)
- plt.gca().invert_xaxis()
- plt.gca().set_aspect("equal")
- plt.legend()
- plt.savefig(f"single_dish_negative_{kk}.png")
- plt.close()
-
- vv = obs.apply_flags(vv.val)
- lim = np.max(np.abs(vv))
- plt.hist(vv.ravel(), bins=np.arange(-lim, lim + 1, 1))
- plt.xlabel("Flux [???]")
- plt.ylabel("Count")
- plt.title("Single dish data")
- plt.savefig(f"single_dish_data_{kk}.png")
- plt.close()
-
-
-def plot_pointings():
- xs, ys = [], []
- for vv in c.interobs.values():
- xs.append(vv.direction.phase_center[0])
- ys.append(vv.direction.phase_center[1])
- xs, ys = np.array(xs), np.array(ys)
- plt.scatter(xs, ys, label="Interferometer")
- for kk, vv in c.sdobs.items():
- plt.scatter(
- *vv.pointings.phase_centers.T, label="Single dish " + kk, alpha=0.2, s=1
- )
- plt.scatter(*c.global_phase_center_casa, label="Phase center CASA")
- plt.scatter(*c.global_phase_center, label="Phase center resolve")
- plt.xlim(
- [c.global_phase_center[0] - c.xfov / 2, c.global_phase_center[0] + c.xfov / 2]
- )
- plt.ylim(
- [c.global_phase_center[1] - c.yfov / 2, c.global_phase_center[1] + c.yfov / 2]
- )
- plt.gca().invert_xaxis()
- plt.gca().set_aspect("equal")
- plt.legend()
- plt.savefig("pointings.png")
- plt.close()
-
-
-def plot_sky_prior():
- p = ift.Plot()
- pspecs = []
- for _ in range(8):
- pos = ift.from_random(c.sky.domain)
- p.add(c.logsky(pos))
- pspecs.append(c.logsky.power_spectrum.force(pos))
- p.add(pspecs)
- p.output(name="sky_prior.png")
-
-
-def main():
- plot_sky_prior()
- plot_pointings()
- plot_sd_data()
-
-
-if __name__ == "__main__":
- main()
diff --git a/misc/convert_alma_data.py b/misc/convert_alma_data.py
deleted file mode 100644
index fd578ae2..00000000
--- a/misc/convert_alma_data.py
+++ /dev/null
@@ -1,141 +0,0 @@
-# SPDX-License-Identifier: GPL-3.0-or-later
-# Copyright(C) 2021 Max-Planck-Society
-# Author: Philipp Arras
-
-from os.path import join, split, splitext
-
-import numpy as np
-from casacore.tables import table, taql
-
-import resolve as rve
-
-
-def cleanup(weight, vis, freq, uvw=None):
- assert weight.ndim == 2
- assert vis.shape == weight.shape
- assert freq.ndim == 1
- assert freq.size == vis.shape[1]
- assert uvw.shape[0] == vis.shape[0]
- # Clean up rows
- ind = np.any(weight != 0.0, axis=1)
- mywgt = weight[ind]
- myvis = vis[ind]
- if uvw is not None:
- myuvw = uvw[ind]
-
- # Clean up freqs
- ind = np.any(mywgt != 0.0, axis=0)
- mywgt = mywgt[:, ind]
- myvis = myvis[:, ind]
- myfreq = freq[ind]
- myvis[mywgt == 0.0] = 0.0
- if uvw is None:
- return mywgt, myvis, myfreq
- return mywgt, myvis, myfreq, myuvw
-
-
-_CASACORE_TABLE_CFG = {"readonly": True, "ack": False}
-minfreq = 230172900000.0
-# if outfbase == "compact":
-# minfreq -= 27146364.99822998
-maxfreq = minfreq * 1.000015
-
-for f, spectral_window, single_dish, fieldid in [
- # ("extended.ms", 0, False, None),
- # ("compact.ms", 0, False, None),
- ("totalpower.ms", 17, True, 1)
-]:
- outfbase = split(splitext(f)[0])[1]
-
- with table(join(f, "SPECTRAL_WINDOW"), **_CASACORE_TABLE_CFG) as t:
- freq = taql("select CHAN_FREQ from $t where ROWID()=$spectral_window").getcol(
- "CHAN_FREQ"
- )
- freq = np.squeeze(freq)
-
- if not np.all(np.diff(freq) >= 0):
- a = freq[::-1]
- assert np.all(np.diff(a) >= 0)
- ind1 = freq.size - np.searchsorted(a, minfreq) + 1
- ind0 = freq.size - np.searchsorted(a, maxfreq) + 1
- else:
- assert np.all(np.diff(freq) >= 0)
- ind0 = np.searchsorted(freq, minfreq)
- ind1 = np.searchsorted(freq, maxfreq)
- channel_slice = slice(ind0, ind1)
- freq = freq[channel_slice]
-
- if not single_dish:
- with table(join(f, "FIELD"), **_CASACORE_TABLE_CFG) as t:
- field_directions = np.squeeze(t.getcol("DELAY_DIR"))
- assert field_directions.ndim == 2
- assert field_directions.shape[1] == 2
-
- with table(f, **_CASACORE_TABLE_CFG) as t:
- sel = lambda x: np.array(
- taql(f"select {x} from $t where DATA_DESC_ID=$spectral_window").getcol(x)
- )
- wgt = sel("WEIGHT")[:, None]
- wgt = np.repeat(wgt, len(freq), axis=1)
- flag = sel("FLAG")[:, channel_slice]
- wgt[flag] = 0.0
- del flag
-
- vis = sel("FLOAT_DATA" if single_dish else "DATA")[:, channel_slice]
- # # FIXME QUESTION Why is there negative data?
- # if single_dish:
- # assert not np.iscomplexobj(vis)
- # ind = vis * np.sqrt(wgt) <= 2.0
- # print(
- # f"{np.sum(ind)} data points are not 2-sigma positive ({np.sum(ind)/ind.size*100:.1f} %)"
- # )
- # wgt[ind] = 0.0
-
- # Average polarization
- assert vis.shape[2] == 2
- assert wgt.shape[2] == 2
- vis = np.sum(wgt * vis, axis=-1)
- wgt = np.sum(wgt, axis=2)
- vis /= wgt
- field = sel("FIELD_ID")
- if single_dish:
- wgt[field != fieldid] = 0.0
- del field
- else:
- uvw = sel("UVW")
-
- if single_dish:
- with table(f, readonly=True) as t:
- pointings = taql(
- "select DIRECTION from $t::POINTING where ROWID() in [select from $t where DATA_DESC_ID=$spectral_window giving [ROWID()]]"
- ).getcol("DIRECTION")
- mywgt, myvis, myfreq, mypointings = cleanup(wgt, vis, freq, pointings)
- assert mypointings.shape[1] == 1
- mypointings = mypointings[:, 0]
- mypointings = rve.Directions(mypointings, 2000)
- obs = rve.SingleDishObservation(
- mypointings,
- np.ascontiguousarray(myvis[None]),
- np.ascontiguousarray(mywgt[None]),
- rve.Polarization.trivial(),
- myfreq,
- )
- obs.save(f"{outfbase}.npz", compress=False)
- else:
- for ifield in set(field):
- ww = wgt.copy()
- ww[field != ifield] = 0.0
- mywgt, myvis, myfreq, myuvw = cleanup(ww, vis, freq, uvw)
- print(ifield, myuvw.shape, myvis.shape, mywgt.shape, myfreq.shape)
-
- obs = rve.Observation(
- rve.AntennaPositions(np.ascontiguousarray(myuvw)),
- np.ascontiguousarray(myvis[None]),
- np.ascontiguousarray(mywgt[None]),
- rve.Polarization.trivial(),
- myfreq,
- rve.Direction(field_directions[ifield], 2000),
- )
- obs.save(f"{outfbase}_field{ifield}.npz", compress=False)
- check = rve.Observation.load(f"{outfbase}_field{ifield}.npz")
- assert check == obs
--
GitLab