Skip to content
Snippets Groups Projects

Mosaicing

Closed Philipp Arras requested to merge parras/resolve:mosaicing into master
Compare and
16 files
+ 921
70
Compare changes
  • Side-by-side
  • Inline
Files
16
+ 113
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 = 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
 
 
# 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():
 
# 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
 
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):
 
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
 
 
 
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
Loading