Mosaicing

demo/cfg_mosaicing.py

+# 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