From e3b6f4b08d1a0d22580c1cc2d04498a68aad9bc1 Mon Sep 17 00:00:00 2001
From: Philipp Arras <c@philipp-arras.de>
Date: Tue, 17 May 2022 13:26:21 +0200
Subject: [PATCH] Work on multires CLEAN
---
playground/multires_clean.py | 14 +++++++++++---
resolve/data/direction.py | 14 ++++++++++++++
resolve/dirty_image.py | 6 +++---
resolve/response.py | 15 ++++++++++-----
resolve/ubik_tools/fits.py | 15 +++++----------
5 files changed, 43 insertions(+), 21 deletions(-)
diff --git a/playground/multires_clean.py b/playground/multires_clean.py
index 5d0e8a26..f5e6f3b3 100644
--- a/playground/multires_clean.py
+++ b/playground/multires_clean.py
@@ -17,6 +17,11 @@ do_wgridding = True
epsilon = 1e-8
nthreads = 6
+if True:
+ do_wgridding = False
+ epsilon = 1e-4
+ nthreads = 6
+
def highres_domain_and_lowrew_mask(dom, oversampling, xmin, xmax, ymin, ymax):
rve.assert_sky_domain(dom)
@@ -44,7 +49,8 @@ def highres_domain_and_lowrew_mask(dom, oversampling, xmin, xmax, ymin, ymax):
mask_lowres[xmin:xmax, ymin:ymax] = 0.
mask_lowres = ift.makeField(sdom, mask_lowres)
- center = ((xmax-xmin)//2 - Nx//2)*Dstx, ((ymax-ymin)//2 - Ny//2)*Dsty
+ center_highres = (xmin + (xmax-xmin)//2 - Nx//2)*Dstx, (ymin + (ymax-ymin)//2 - Ny//2)*Dsty
+ print(np.array(center_highres) / rve.ARCMIN2RAD)
npix = np.array([(xmax-xmin)*oversampling, (ymax-ymin)*oversampling])
assert (npix % 2 == 0).all()
@@ -75,6 +81,8 @@ def main():
squeeze = ift.ContractionOperator(dom, (0, 1, 2))
dom_highres, mask_lowres, center_highres = highres_domain_and_lowrew_mask(dom, oversampling_highres, "-0.94arcmin", "-0.78arcmin", "-0.40arcmin", "-0.28arcmin")
+ dom_highres, mask_lowres, center_highres = highres_domain_and_lowrew_mask(dom, oversampling_highres, "-0.2arcmin", "0.2arcmin", "-0.2arcmin", "0.2arcmin")
+ direction_highres = rve.Direction(center_highres, ms.direction.equinox) + ms.direction
R = rve.InterferometryResponse(ms, dom, do_wgridding, epsilon, nthreads=nthreads)
R_highres = rve.InterferometryResponse(ms, dom_highres, do_wgridding, epsilon,
@@ -85,8 +93,8 @@ def main():
psf = get_psf(ms, dom)
psf_highres = get_psf(ms, dom, center=center_highres)
- rve.ubik_tools.field2fits(dirty, "dirty.fits")
- rve.ubik_tools.field2fits(dirty_highres, "dirty_highres.fits")
+ rve.ubik_tools.field2fits(dirty, "dirty.fits", ms.direction)
+ rve.ubik_tools.field2fits(dirty_highres, "dirty_highres.fits", direction_highres)
exit()
rve.ubik_tools.field2fits(psf, "psf.fits")
diff --git a/resolve/data/direction.py b/resolve/data/direction.py
index f17ea2c9..6faf0d3d 100644
--- a/resolve/data/direction.py
+++ b/resolve/data/direction.py
@@ -12,8 +12,11 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2020-2021 Max-Planck-Society
+# Copyright(C) 2022 Max-Planck-Society, Philipp Arras
# Author: Philipp Arras
+import numpy as np
+
from ..util import compare_attributes, my_asserteq
@@ -55,6 +58,17 @@ class Direction:
return False
return compare_attributes(self, other, ("_pc", "_e"))
+ def __add__(self, other):
+ if self.equinox != other.equinox:
+ s = "Equinox need to be the same. Got {self.equinox} and {other.equinox}"
+ raise ValueError(s)
+ pc0, pc1 = np.array(self.phase_center), np.array(other.phase_center)
+ return Direction(tuple(pc0 + pc1), self.equinox)
+
+ def __sub__(self, other):
+ pc = np.array((self.phase_center[0], self.phase_center[1]))
+ return Direction(-pc, self.equinox) + other
+
class Directions:
def __init__(self, phase_centers, equinox):
diff --git a/resolve/dirty_image.py b/resolve/dirty_image.py
index 76337c27..3cbfc3e0 100644
--- a/resolve/dirty_image.py
+++ b/resolve/dirty_image.py
@@ -26,10 +26,10 @@ from .util import assert_sky_domain
def dirty_image(observation, weighting, sky_domain, do_wgridding, epsilon,
- vis=None, weight=None, nthreads=1):
+ vis=None, weight=None, center=(0., 0.), nthreads=1):
assert_sky_domain(sky_domain)
R = InterferometryResponse(observation, sky_domain, do_wgridding=do_wgridding,
- epsilon=epsilon, nthreads=nthreads)
+ epsilon=epsilon, nthreads=nthreads, center=center)
w = observation.weight if weight is None else weight
d = observation.vis if vis is None else vis
vol = sky_domain[-1].scalar_dvol
@@ -37,7 +37,7 @@ def dirty_image(observation, weighting, sky_domain, do_wgridding, epsilon,
return R.adjoint(d * w/w.s_sum()) / vol**2
elif weighting == "uniform":
w = uniform_weights(observation, sky_domain)
- return R.adjoint(d * w/w.s_sum()) / vol**2
+ return R.adjoint(d * w/w.s_sum()) / vol**2 * np.sqrt(vol)**2
raise RuntimeError
diff --git a/resolve/response.py b/resolve/response.py
index 4103049a..5baa8b95 100644
--- a/resolve/response.py
+++ b/resolve/response.py
@@ -29,8 +29,10 @@ from .dtype_converter import DtypeConverter
from .logger import logger
-def InterferometryResponse(observation, domain, do_wgridding, epsilon, verbosity=0, nthreads=1):
- R = _InterferometryResponse(observation, domain, do_wgridding, epsilon, verbosity, nthreads=nthreads)
+def InterferometryResponse(observation, domain, do_wgridding, epsilon, verbosity=0, nthreads=1,
+ center=(0., 0.)):
+ R = _InterferometryResponse(observation, domain, do_wgridding, epsilon, verbosity, center,
+ nthreads=nthreads)
pol = polarization_converter(R.target, observation.vis.domain)
if observation.double_precision:
dtype = DtypeConverter(domain, np.float64, np.float64, "Response input") # do nothing
@@ -45,7 +47,7 @@ def InterferometryResponse(observation, domain, do_wgridding, epsilon, verbosity
class _InterferometryResponse(ift.LinearOperator):
- def __init__(self, observation, domain, do_wgridding, epsilon, verbosity, nthreads=1):
+ def __init__(self, observation, domain, do_wgridding, epsilon, verbosity, center, nthreads=1):
assert isinstance(observation, Observation)
domain = ift.DomainTuple.make(domain)
assert_sky_domain(domain)
@@ -86,7 +88,8 @@ class _InterferometryResponse(ift.LinearOperator):
# FIXME Include mask in the future here? Problem is that
# the mask may be different for different polarizations
rrr = SingleResponse(domain[3], ooo.uvw, ooo.freq, do_wgridding=do_wgridding,
- epsilon=epsilon, verbosity=verbosity, nthreads=nthreads)
+ epsilon=epsilon, verbosity=verbosity, nthreads=nthreads,
+ center=center)
sr_tmp.append(rrr)
t_tmp.append(tind)
f_tmp.append(find)
@@ -138,7 +141,7 @@ class _InterferometryResponse(ift.LinearOperator):
class SingleResponse(ift.LinearOperator):
def __init__(self, domain, uvw, freq, do_wgridding, epsilon, mask=None, facets=(1, 1),
- verbosity=0, nthreads=1):
+ center=(0., 0.), verbosity=0, nthreads=1):
my_assert_isinstance(facets, tuple)
my_assert_isinstance(do_wgridding, bool)
my_assert_isinstance(epsilon, float)
@@ -160,6 +163,8 @@ class SingleResponse(ift.LinearOperator):
"do_wgridding": do_wgridding,
"nthreads": nthreads,
"flip_v": True,
+ "center_x": center[0],
+ "center_y": center[1],
"verbosity": verbosity
}
if mask is not None:
diff --git a/resolve/ubik_tools/fits.py b/resolve/ubik_tools/fits.py
index eb8ee5e8..3fab3b28 100644
--- a/resolve/ubik_tools/fits.py
+++ b/resolve/ubik_tools/fits.py
@@ -22,21 +22,16 @@ import numpy as np
from ..constants import DEG2RAD
from ..util import assert_sky_domain
-from ..data.observation import Observation
+from ..data.direction import Direction
-def field2fits(field, file_name, observations=[]):
+def field2fits(field, file_name, direction=None):
import astropy.io.fits as pyfits
from astropy.time import Time
- if isinstance(observations, Observation):
- observations = [observations]
- if len(observations) == 0:
- direction = None
- else:
- if len(observations) > 1:
- warn("field2fits: Include only info of first observation into fits file.")
- direction = observations[0].direction
+ if not isinstance(direction, Direction):
+ s = f"`direction` needs to be an instance of `rve.Direction`. Got {direction}"
+ raise TypeError(s)
dom = field.domain
assert_sky_domain(dom)
--
GitLab