From 0174bf27ab9b4b4cd33530b78b4ddc616ecbcc74 Mon Sep 17 00:00:00 2001 From: paolomassa Date: Tue, 28 May 2024 09:39:51 +0200 Subject: [PATCH 01/12] Fixed bug in the default creation of VisMeta --- xrayvision/visibility.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/xrayvision/visibility.py b/xrayvision/visibility.py index 0740c07..f7e4dfe 100644 --- a/xrayvision/visibility.py +++ b/xrayvision/visibility.py @@ -285,10 +285,6 @@ def __init__( self._uv_key = "uv" self._units_key = "units" - # Build meta. Make sure that phase center is included. - if not isinstance(meta, VisMetaABC): - meta = VisMeta(meta) - # Construct underlying data object. dims = [f"dim{i}" for i in range(0, len(visibilities.shape))] dims[_uv_axis] = self._uv_key @@ -307,7 +303,7 @@ def __init__( if phase_uncertainty is not None: data[self._phase_uncert_key] = (dims, phase_uncertainty.to_value(phase.unit)) if meta is None: - meta = VisMeta(dict()) + meta = VisMeta() vis_labels = getattr(meta, "vis_labels", None) if vis_labels is not None: if len(vis_labels) != nvis: From ac782c71eb9e2d6ed6826085f1acb743c7def2c6 Mon Sep 17 00:00:00 2001 From: paolomassa Date: Tue, 28 May 2024 10:05:56 +0200 Subject: [PATCH 02/12] Add changelog file --- changelog/66.bugfix.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 changelog/66.bugfix.rst diff --git a/changelog/66.bugfix.rst b/changelog/66.bugfix.rst new file mode 100644 index 0000000..3f23915 --- /dev/null +++ b/changelog/66.bugfix.rst @@ -0,0 +1 @@ +Fix bug that happens when you create a Visibility object with default meta. From 243ace1c2353987518d35cb38917fc7ecf7de7a9 Mon Sep 17 00:00:00 2001 From: paolomassa Date: Wed, 29 May 2024 12:43:59 +0200 Subject: [PATCH 03/12] Made the software compatible with the Visibilities class --- examples/stix.py | 2 +- xrayvision/clean.py | 8 +- xrayvision/imaging.py | 50 ++++++------ xrayvision/mem.py | 22 +++--- xrayvision/tests/test_imaging.py | 41 +++++----- xrayvision/tests/test_mem.py | 2 +- xrayvision/tests/test_transform.py | 114 ++++++++++++++-------------- xrayvision/tests/test_visibility.py | 14 +++- xrayvision/transform.py | 42 +++++----- xrayvision/visibility.py | 2 +- 10 files changed, 154 insertions(+), 143 deletions(-) diff --git a/examples/stix.py b/examples/stix.py index c8511c3..5b908e6 100644 --- a/examples/stix.py +++ b/examples/stix.py @@ -26,7 +26,7 @@ stix_data = pickle.load(urllib.request.urlopen("https://pub099.cs.technik.fhnw.ch/demo/stix_vis.pkl")) time_range, energy_range, offset, stix_vis = stix_data -stix_vis.phase_centre = [0, 0] * apu.arcsec +stix_vis.phase_center = [0, 0] * apu.arcsec stix_vis.offset = offset ############################################################################### diff --git a/xrayvision/clean.py b/xrayvision/clean.py index 32327a3..e6f8f9b 100644 --- a/xrayvision/clean.py +++ b/xrayvision/clean.py @@ -21,7 +21,7 @@ from sunpy.map.map_factory import Map from xrayvision.imaging import vis_psf_image, vis_to_map -from xrayvision.visibility import Visibility +from xrayvision.visibility import Visibilities __all__ = ["clean", "vis_clean", "ms_clean", "vis_ms_clean"] @@ -98,7 +98,7 @@ def clean( # raise ValueError('') pad = [0 if x % 2 == 0 else 1 for x in dirty_map.shape] - # Assume beam, map phase_centre is in middle + # Assume beam, map phase_center is in middle beam_center = (dirty_beam.shape[0] - 1) / 2.0, (dirty_beam.shape[1] - 1) / 2.0 map_center = (dirty_map.shape[0] - 1) / 2.0, (dirty_map.shape[1] - 1) / 2.0 @@ -173,7 +173,7 @@ def clean( @u.quantity_input def vis_clean( - vis: Visibility, + vis: Visibilities, shape: Quantity[u.pix], pixel_size: Quantity[u.arcsec / u.pix], clean_beam_width: Optional[Quantity[u.arcsec]] = 4.0, @@ -403,7 +403,7 @@ def ms_clean( def vis_ms_clean( - vis: Visibility, + vis: Visibilities, shape: Quantity[u.pix], pixel_size: Quantity[u.arcsec / u.pix], scales: Optional[Iterable], diff --git a/xrayvision/imaging.py b/xrayvision/imaging.py index c304266..d23de65 100644 --- a/xrayvision/imaging.py +++ b/xrayvision/imaging.py @@ -6,7 +6,7 @@ from sunpy.map import GenericMap, Map from xrayvision.transform import dft_map, idft_map -from xrayvision.visibility import Visibility +from xrayvision.visibility import Visibility, Visibilities __all__ = [ "get_weights", @@ -24,7 +24,7 @@ WEIGHT_SCHEMES = ("natural", "uniform") -def get_weights(vis: Visibility, scheme: str = "natural", norm: bool = True) -> np.ndarray: +def get_weights(vis: Visibilities, scheme: str = "natural", norm: bool = True) -> np.ndarray: r""" Return spatial frequency weight factors for each visibility. @@ -47,7 +47,7 @@ def get_weights(vis: Visibility, scheme: str = "natural", norm: bool = True) -> raise ValueError(f"Invalid weighting scheme {scheme}, must be one of: {WEIGHT_SCHEMES}") weights = np.sqrt(vis.u**2 + vis.v**2).value if scheme == "natural": - weights = np.ones_like(vis.vis, dtype=float) + weights = np.ones_like(vis.visibilities, dtype=float) if norm: weights /= weights.sum() @@ -97,11 +97,11 @@ def image_to_vis( *, u: Quantity[apu.arcsec**-1], v: Quantity[apu.arcsec**-1], - phase_centre: Optional[Quantity[apu.arcsec]] = (0.0, 0.0) * apu.arcsec, + phase_center: Optional[Quantity[apu.arcsec]] = (0.0, 0.0) * apu.arcsec, pixel_size: Optional[Quantity[apu.arcsec / apu.pix]] = 1.0 * apu.arcsec / apu.pix, -) -> Visibility: +) -> Visibilities: r""" - Return a Visibility created from the image and u, v sampling. + Return a Visibilities object created from the image and u, v sampling. Parameters ---------- @@ -111,27 +111,27 @@ def image_to_vis( Array of u coordinates where the visibilities will be evaluated v : Array of v coordinates where the visibilities will be evaluated - phase_centre : - The coordinates the phase_centre. + phase_center : + The coordinates the phase_center. pixel_size : Size of pixels, if only one value is passed, assume square pixels (repeating the value). Returns ------- : - The new visibility object + The new Visibilities object """ pixel_size = validate_and_expand_kwarg(pixel_size, "pixel_size") if not (apu.get_physical_type((1 / u).unit) == ANGLE and apu.get_physical_type((1 / v).unit) == ANGLE): raise ValueError("u and v must be inverse angle (e.g. 1/deg or 1/arcsec") - vis = dft_map(image, u=u, v=v, phase_centre=phase_centre, pixel_size=pixel_size) - return Visibility(vis, u=u, v=v, offset=phase_centre) + vis = dft_map(image, u=u, v=v, phase_center=phase_center, pixel_size=pixel_size) + return Visibilities(vis, u=u, v=v, phase_center=phase_center) @apu.quantity_input() def vis_to_image( - vis: Visibility, + vis: Visibilities, shape: Quantity[apu.pix] = (65, 65) * apu.pixel, pixel_size: Optional[Quantity[apu.arcsec / apu.pix]] = 1 * apu.arcsec / apu.pix, scheme: str = "natural", @@ -161,7 +161,7 @@ def vis_to_image( shape = shape.to(apu.pixel) weights = get_weights(vis, scheme=scheme) bp_arr = idft_map( - vis.vis, u=vis.u, v=vis.v, shape=shape, weights=weights, pixel_size=pixel_size, phase_centre=vis.phase_centre + vis.visibilities, u=vis.u, v=vis.v, shape=shape, weights=weights, pixel_size=pixel_size, phase_center=vis.phase_center ) return bp_arr @@ -169,7 +169,7 @@ def vis_to_image( @apu.quantity_input def vis_psf_map( - vis: Visibility, + vis: Visibilities, *, shape: Quantity[apu.pix] = (65, 65) * apu.pixel, pixel_size: Optional[Quantity[apu.arcsec / apu.pix]] = 1 * apu.arcsec / apu.pix, @@ -203,7 +203,7 @@ def vis_psf_map( @apu.quantity_input() def vis_psf_image( - vis: Visibility, + vis: Visibilities, *, shape: Quantity[apu.pix] = (65, 65) * apu.pixel, pixel_size: Quantity[apu.arcsec / apu.pix] = 1 * apu.arcsec / apu.pix, @@ -237,14 +237,14 @@ def vis_psf_image( # Make sure psf is always odd so power is in exactly one pixel shape = [s // 2 * 2 + 1 for s in shape.to_value(apu.pix)] * shape.unit psf_arr = idft_map( - np.ones(vis.vis.shape) * vis.vis.unit, u=vis.u, v=vis.v, shape=shape, weights=weights, pixel_size=pixel_size + np.ones(vis.visibilities.shape) * vis.visibilities.unit, u=vis.u, v=vis.v, shape=shape, weights=weights, pixel_size=pixel_size ) return psf_arr @apu.quantity_input() def vis_to_map( - vis: Visibility, + vis: Visibilities, shape: Quantity[apu.pix] = (65, 65) * apu.pixel, pixel_size: Optional[Quantity[apu.arcsec / apu.pix]] = 1 * apu.arcsec / apu.pixel, scheme: Optional[str] = "natural", @@ -278,7 +278,7 @@ def vis_to_map( @apu.quantity_input() -def generate_header(vis: Visibility, *, shape: Quantity[apu.pix], pixel_size: Quantity[apu.arcsec / apu.pix]) -> dict: +def generate_header(vis: Visibilities, *, shape: Quantity[apu.pix], pixel_size: Quantity[apu.arcsec / apu.pix]) -> dict: r""" Generate a map head given the visibilities, pixel size and shape @@ -296,8 +296,8 @@ def generate_header(vis: Visibility, *, shape: Quantity[apu.pix], pixel_size: Qu : """ header = { - "crval1": (vis.offset[1]).to_value(apu.arcsec), - "crval2": (vis.offset[0]).to_value(apu.arcsec), + "crval1": (vis.phase_center[1]).to_value(apu.arcsec), + "crval2": (vis.phase_center[0]).to_value(apu.arcsec), "cdelt1": (pixel_size[1] * apu.pix).to_value(apu.arcsec), "cdelt2": (pixel_size[0] * apu.pix).to_value(apu.arcsec), "ctype1": "HPLN-TAN", @@ -312,9 +312,9 @@ def generate_header(vis: Visibility, *, shape: Quantity[apu.pix], pixel_size: Qu @apu.quantity_input() -def map_to_vis(amap: GenericMap, *, u: Quantity[1 / apu.arcsec], v: Quantity[1 / apu.arcsec]) -> Visibility: +def map_to_vis(amap: GenericMap, *, u: Quantity[1 / apu.arcsec], v: Quantity[1 / apu.arcsec]) -> Visibilities: r""" - Return a Visibility object created from the map, sampling it at give `u`, `v` coordinates. + Return a Visibilities object created from the map, sampling it at give `u`, `v` coordinates. Parameters ---------- @@ -328,7 +328,7 @@ def map_to_vis(amap: GenericMap, *, u: Quantity[1 / apu.arcsec], v: Quantity[1 / Returns ------- : - The new visibility object + The new Visibilities object """ if not apu.get_physical_type(1 / u) == ANGLE and apu.get_physical_type(1 / v) == ANGLE: @@ -346,6 +346,6 @@ def map_to_vis(amap: GenericMap, *, u: Quantity[1 / apu.arcsec], v: Quantity[1 / if "cdelt2" in meta: new_psize[0] = float(meta["cdelt2"]) - vis = image_to_vis(amap.data, u=u, v=v, pixel_size=new_psize * apu.arcsec / apu.pix) - vis.offset = new_pos * apu.arcsec + vis = image_to_vis(amap.quantity, u=u, v=v, pixel_size=new_psize * apu.arcsec / apu.pix, phase_center=new_pos * apu.arcsec) + return vis diff --git a/xrayvision/mem.py b/xrayvision/mem.py index 485e86f..3c68470 100644 --- a/xrayvision/mem.py +++ b/xrayvision/mem.py @@ -27,7 +27,7 @@ "mem", ] -from xrayvision.visibility import Visibility +from xrayvision.visibility import Visibility, Visibilities logger = get_logger(__name__, "DEBUG") @@ -81,8 +81,8 @@ def _get_fourier_matrix(vis, shape=(64, 64) * apu.pix, pixel_size=(4.0312500, 4. The complex Fourier matrix """ m, n = shape.to("pix") - y = generate_xy(m, phase_centre=0 * apu.arcsec, pixel_size=pixel_size[1]) - x = generate_xy(n, phase_centre=0 * apu.arcsec, pixel_size=pixel_size[0]) + y = generate_xy(m, phase_center=0 * apu.arcsec, pixel_size=pixel_size[1]) + x = generate_xy(n, phase_center=0 * apu.arcsec, pixel_size=pixel_size[0]) x, y = np.meshgrid(x, y, indexing="ij") uv = np.vstack([vis.u, vis.v]) # Check apu are correct for exp need to be dimensionless and then remove apu for speed @@ -192,8 +192,8 @@ def _prepare_for_optimise(pixel, shape, vis): # 'Hv' is composed by the union of its real and imaginary part Hv = np.concatenate([ReHv, ImHv], axis=-1) # Division of real and imaginary part of the visibilities - ReV = np.real(vis.vis) - ImV = np.imag(vis.vis) + ReV = np.real(vis.visibilities) + ImV = np.imag(vis.visibilities) # 'Visib' is composed by the real and imaginary part of the visibilities Visib = np.concatenate([ReV, ImV], axis=-1) # Standard deviation of the real and imaginary part of the visibilities @@ -256,7 +256,7 @@ def _get_mean_visibilities(vis, shape, pixel): v = np.zeros(n_vis) * (1 / apu.arcsec) den = np.zeros(n_vis) weights = np.zeros_like(vis.amplitude_error**2) - visib = np.zeros_like(vis.vis) + visib = np.zeros_like(vis.visibilities) for ip in range(n_vis): # what about 0.5 pix offset i = ru[ip].to_value("pix").astype(int) @@ -269,7 +269,7 @@ def _get_mean_visibilities(vis, shape, pixel): # we save in 'u' and 'v' the u and v coordinates of the first frequency that corresponds # to the position (i, j) of the discretization of the (u,v)-plane 'iuarr' - visib[count] = vis.vis[ip] + visib[count] = vis.visibilities[ip] weights[count] = vis.amplitude_error[ip] ** 2.0 den[count] = 1.0 iuarr[i, j] = count @@ -277,7 +277,7 @@ def _get_mean_visibilities(vis, shape, pixel): count += 1 else: # Save the sum of the visibilities that correspond to the same position (i, j) - visib[iuarr[i, j].astype(int)] += vis.vis[ip] + visib[iuarr[i, j].astype(int)] += vis.visibilities[ip] # Save the number of the visibilities that correspond to the same position (i, j) den[iuarr[i, j].astype(int)] += 1.0 # Save the sum of the variances of the amplitudes of the visibilities that @@ -296,7 +296,7 @@ def _get_mean_visibilities(vis, shape, pixel): # correspond to the same position in the discretization of the (u,v)-plane weights = np.sqrt(weights[:count]) / den - return SimpleNamespace(u=u, v=v, vis=visib, amplitude_error=weights) + return SimpleNamespace(u=u, v=v, visibilities=visib, amplitude_error=weights) def _proximal_entropy(y, m, lamba, Lip, tol=10**-10): @@ -601,7 +601,7 @@ def _get_percent_lambda(vis): if ibig.size < 2: snr_value = -1 else: - snr_value, _ = resistant_mean((np.abs(vis.vis[ibig]) / vis.amplitude_error[ibig]).flatten(), 3) + snr_value, _ = resistant_mean((np.abs(vis.visibilities[ibig]) / vis.amplitude_error[ibig]).flatten(), 3) # TODO magic numbers percent_lambda = 2 / (snr_value**2 + 90) @@ -611,7 +611,7 @@ def _get_percent_lambda(vis): @apu.quantity_input def mem( - vis: Visibility, + vis: Visibilities, shape: Quantity[apu.pix], pixel_size: Quantity[apu.arcsec / apu.pix], *, diff --git a/xrayvision/tests/test_imaging.py b/xrayvision/tests/test_imaging.py index 82819c4..d19981f 100644 --- a/xrayvision/tests/test_imaging.py +++ b/xrayvision/tests/test_imaging.py @@ -7,7 +7,7 @@ from xrayvision.imaging import image_to_vis, map_to_vis, vis_psf_image, vis_to_image, vis_to_map from xrayvision.transform import dft_map, generate_uv, idft_map -from xrayvision.visibility import Visibility +from xrayvision.visibility import Visibility, Visibilities @pytest.fixture @@ -93,7 +93,7 @@ def test_vis_to_psf(pixel_size, uv): weights = np.sqrt(u**2 + v**2).value weights /= weights.sum() psf_calc = idft_map(obs_vis, u=u, v=v, shape=[65, 65] * apu.pix, pixel_size=ps, weights=weights) - vis = Visibility(obs_vis, u=u, v=v) + vis = Visibilities(obs_vis, u, v) res = vis_psf_image(vis, shape=[65, 65] * apu.pixel, pixel_size=ps, scheme="uniform") assert_allclose(res, psf_calc) @@ -108,7 +108,7 @@ def test_vis_to_image_against_idft(uv): bp_calc = idft_map( obs_vis, u=u, v=v, shape=[65, 65] * apu.pix, pixel_size=[2, 2] * apu.arcsec / apu.pix, weights=weights ) - vis = Visibility(obs_vis, u=u, v=v) + vis = Visibilities(obs_vis, u, v) res = vis_to_image(vis, shape=[65, 65] * apu.pixel, pixel_size=2 * apu.arcsec / apu.pix, scheme="uniform") assert np.allclose(bp_calc, res) @@ -117,7 +117,7 @@ def test_image_to_vis(): m = n = 33 size = m * n # Set up empty map - image = np.zeros((m, n)) + image = np.zeros((m, n)) * apu.ct / apu.arcsec**2 # Calculate full u, v coverage so will be equivalent to a discrete Fourier transform (DFT) u = generate_uv(m * apu.pix) @@ -127,15 +127,15 @@ def test_image_to_vis(): # For an empty map visibilities should all be zero (0+0j) empty_vis = image_to_vis(image, u=v, v=v) - assert np.array_equal(empty_vis.phase_centre, (0.0, 0.0) * apu.arcsec) - assert np.array_equal(empty_vis.vis, np.zeros(n * m, dtype=complex)) + assert np.array_equal(empty_vis.phase_center, (0.0, 0.0) * apu.arcsec) + assert np.array_equal(empty_vis.visibilities, np.zeros(n * m, dtype=complex)) def test_image_to_vis_center(): m = n = 33 size = m * n # Set up empty map - image = np.zeros((m, n)) + image = np.zeros((m, n)) * apu.ct / apu.arcsec**2 # Calculate full u, v coverage so will be equivalent to a discrete Fourier transform (DFT) u = generate_uv(m * apu.pix) @@ -144,9 +144,9 @@ def test_image_to_vis_center(): u, v = np.array([u, v]).reshape(2, size) / apu.arcsec # For an empty map visibilities should all be zero (0+0j) - empty_vis = image_to_vis(image, u=u, v=v, phase_centre=(2.0, -3.0) * apu.arcsec) - assert np.array_equal(empty_vis.offset, (2.0, -3.0) * apu.arcsec) - assert np.array_equal(empty_vis.vis, np.zeros(n * m, dtype=complex)) + empty_vis = image_to_vis(image, u=u, v=v, phase_center=(2.0, -3.0) * apu.arcsec) + assert np.array_equal(empty_vis.phase_center, (2.0, -3.0) * apu.arcsec) + assert np.array_equal(empty_vis.visibilities, np.zeros(n * m, dtype=complex)) @pytest.mark.parametrize( @@ -161,8 +161,8 @@ def test_map_to_vis(pos, pixel): pixel = pixel * apu.arcsec / apu.pix # Calculate full u, v coverage so will be equivalent to a discrete Fourier transform (DFT) - u = generate_uv(m * apu.pix, phase_centre=pos[0], pixel_size=pixel[0]) - v = generate_uv(n * apu.pix, phase_centre=pos[1], pixel_size=pixel[1]) + u = generate_uv(m * apu.pix, phase_center=pos[0], pixel_size=pixel[0]) + v = generate_uv(n * apu.pix, phase_center=pos[1], pixel_size=pixel[1]) u, v = np.meshgrid(u, v, indexing="ij") u, v = np.array([u, v]).reshape(2, size) / apu.arcsec @@ -173,14 +173,15 @@ def test_map_to_vis(pos, pixel): "cdelt2": pixel[0].value, "cunit1": "arcsec", "cunit2": "arcsec", + "bunit": "ct / arcsec^2", } # Astropy index order is opposite to that of numpy, is 1st dim is across second down - data = Gaussian2DKernel(6, x_size=n, y_size=m).array + data = Gaussian2DKernel(6, x_size=n, y_size=m).array * apu.ct / apu.arcsec**2 mp = Map((data, header)) vis = map_to_vis(mp, u=u, v=v) - assert np.array_equal(vis.offset, pos) + assert np.array_equal(vis.phase_center, pos) res = vis_to_image(vis, shape=(m, n) * apu.pixel, pixel_size=pixel) assert np.allclose(res, data) @@ -198,9 +199,7 @@ def test_vis_to_image(): u, v = uv # Astropy index order is opposite to that of numpy, is 1st dim is across second down - data = Gaussian2DKernel(6, x_size=n, y_size=m).array - data = np.zeros((m, n)) - data[m // 2, n // 2] = 1 + data = Gaussian2DKernel(6, x_size=n, y_size=m).array * apu.ct / apu.arcsec**2 vis = image_to_vis(data, u=u, v=v) res = vis_to_image(vis, shape=(m, n) * apu.pixel) @@ -219,7 +218,7 @@ def test_vis_to_image_single_pixel_size(): uv = np.array([u, v]).reshape(2, size) / apu.arcsec u, v = uv # Astropy index order is opposite to that of numpy, is 1st dim is across second down - data = Gaussian2DKernel(6, x_size=n, y_size=m).array + data = Gaussian2DKernel(6, x_size=n, y_size=m).array * apu.ct / apu.arcsec**2 vis = image_to_vis(data, u=u, v=v, pixel_size=(2.0, 2.0) * apu.arcsec / apu.pix) res = vis_to_image(vis, shape=(m, n) * apu.pixel, pixel_size=2.0 * apu.arcsec / apu.pix) @@ -239,7 +238,7 @@ def test_vis_to_image_invalid_pixel_size(): u, v = uv # Astropy index order is opposite to that of numpy, is 1st dim is across second down - data = Gaussian2DKernel(6, x_size=n, y_size=m).array + data = Gaussian2DKernel(6, x_size=n, y_size=m).array * apu.ct / apu.arcsec**2 vis = image_to_vis(data, u=u, v=v) with pytest.raises(ValueError): @@ -256,8 +255,8 @@ def test_vis_to_image_invalid_pixel_size(): def test_vis_to_map(m, n, pos, pixel): pos = pos * apu.arcsec pixel = pixel * apu.arcsec / apu.pix - u = generate_uv(m * apu.pix, phase_centre=pos[0], pixel_size=pixel[0]) - v = generate_uv(n * apu.pix, phase_centre=pos[1], pixel_size=pixel[1]) + u = generate_uv(m * apu.pix, phase_center=pos[0], pixel_size=pixel[0]) + v = generate_uv(n * apu.pix, phase_center=pos[1], pixel_size=pixel[1]) u, v = np.meshgrid(u, v, indexing="ij") uv = np.array([u, v]).reshape(2, m * n) / apu.arcsec u, v = uv diff --git a/xrayvision/tests/test_mem.py b/xrayvision/tests/test_mem.py index f2568a5..69a41ee 100644 --- a/xrayvision/tests/test_mem.py +++ b/xrayvision/tests/test_mem.py @@ -57,7 +57,7 @@ def test_mem(): # sub_uv = np.hstack([sub_uv, np.zeros((2, 1))]) / u.arcsec vis = image_to_vis(data * u.dimensionless_unscaled, u=sub_uv[0, :], v=sub_uv[1, :], pixel_size=2 * u.arcsec / u.pix) - setattr(vis, "amplitude_error", np.sqrt(np.abs(vis.vis))) + setattr(vis, "amplitude_error", np.sqrt(np.abs(vis.visibilities))) setattr(vis, "label", [str(x) for x in np.sqrt(x**2 + y**2).flatten()]) res = mem( diff --git a/xrayvision/tests/test_transform.py b/xrayvision/tests/test_transform.py index 4fcb6d8..b08032f 100644 --- a/xrayvision/tests/test_transform.py +++ b/xrayvision/tests/test_transform.py @@ -18,31 +18,31 @@ def test_generate_xy_pixel_size(pixel_size): assert np.array_equal(even, generate_xy(32 * apu.pix, pixel_size=pixel_size)) -@pytest.mark.parametrize("phase_centre", [0, +5.5, -5.5]) -def test_generate_xy_offset(phase_centre): - phase_centre = phase_centre * apu.arcsec - even = np.linspace(-15.5, 15.5, 32) * apu.arcsec + phase_centre - odd = np.linspace(-16, 16, 33) * apu.arcsec + phase_centre +@pytest.mark.parametrize("phase_center", [0, +5.5, -5.5]) +def test_generate_xy_offset(phase_center): + phase_center = phase_center * apu.arcsec + even = np.linspace(-15.5, 15.5, 32) * apu.arcsec + phase_center + odd = np.linspace(-16, 16, 33) * apu.arcsec + phase_center - assert np.array_equal(even, generate_xy(32 * apu.pix, phase_centre=phase_centre)) - assert np.array_equal(odd, generate_xy(33 * apu.pix, phase_centre=phase_centre)) + assert np.array_equal(even, generate_xy(32 * apu.pix, phase_center=phase_center)) + assert np.array_equal(odd, generate_xy(33 * apu.pix, phase_center=phase_center)) @pytest.mark.parametrize( - "phase_centre, pixel_size", [(0, (0.5, 1, 2, 3)), (+5.5, (0.5, 1, 2, 3)), (-5.5, (0.5, 1, 2, 3))] + "phase_center, pixel_size", [(0, (0.5, 1, 2, 3)), (+5.5, (0.5, 1, 2, 3)), (-5.5, (0.5, 1, 2, 3))] ) -def test_generate_xy_offset_size(phase_centre, pixel_size): - phase_centre = phase_centre * apu.arcsec +def test_generate_xy_offset_size(phase_center, pixel_size): + phase_center = phase_center * apu.arcsec pixel_size = pixel_size * apu.arcsec / apu.pix # No sure if this is a good idea could be hard to debug for size in pixel_size: # Odd - odd = np.linspace(-16 * size, 16 * size, 33) * apu.pix + phase_centre - assert np.array_equal(odd, generate_xy(33 * apu.pix, phase_centre=phase_centre, pixel_size=size)) + odd = np.linspace(-16 * size, 16 * size, 33) * apu.pix + phase_center + assert np.array_equal(odd, generate_xy(33 * apu.pix, phase_center=phase_center, pixel_size=size)) # Even - even = np.linspace(-15.5 * size, 15.5 * size, 32) * apu.pix + phase_centre - assert np.array_equal(even, generate_xy(32 * apu.pix, phase_centre=phase_centre, pixel_size=size)) + even = np.linspace(-15.5 * size, 15.5 * size, 32) * apu.pix + phase_center + assert np.array_equal(even, generate_xy(32 * apu.pix, phase_center=phase_center, pixel_size=size)) @pytest.mark.parametrize("pixel_size", [0.5, 1, 2, 3]) @@ -62,30 +62,30 @@ def test_generate_uv_pixel_size(pixel_size): assert np.allclose(even, generate_uv(32 * apu.pix, pixel_size=pixel_size), atol=1e-8 / apu.arcsec) -@pytest.mark.parametrize("phase_centre", [0.0, -5.5, 5.5]) -def test_generate_uv_pixel_offset(phase_centre): - phase_centre = phase_centre * apu.arcsec +@pytest.mark.parametrize("phase_center", [0.0, -5.5, 5.5]) +def test_generate_uv_pixel_offset(phase_center): + phase_center = phase_center * apu.arcsec m = 33 n = 32 # Odd odd = np.linspace(-((m - 1) / 2) * (1 / m), ((m - 1) / 2) * (1 / m), m) / apu.arcsec - if phase_centre.value != 0: - odd += 1 / phase_centre - assert np.allclose(odd, generate_uv(m * apu.pix, phase_centre=phase_centre), atol=1e-8 / apu.arcsec) + if phase_center.value != 0: + odd += 1 / phase_center + assert np.allclose(odd, generate_uv(m * apu.pix, phase_center=phase_center), atol=1e-8 / apu.arcsec) # Even even = (np.arange(n) - n / 2 + 0.5) * (1 / n) / apu.arcsec - if phase_centre.value != 0: - even += 1 / phase_centre - assert np.allclose(even, generate_uv(n * apu.pix, phase_centre=phase_centre), atol=1e-8 / apu.arcsec) + if phase_center.value != 0: + even += 1 / phase_center + assert np.allclose(even, generate_uv(n * apu.pix, phase_center=phase_center), atol=1e-8 / apu.arcsec) @pytest.mark.parametrize( - "phase_centre, pixel_size", [(0, (0.5, 1, 2, 3)), (+5.5, (0.5, 1, 2, 3)), (-5.5, (0.5, 1, 2, 3))] + "phase_center, pixel_size", [(0, (0.5, 1, 2, 3)), (+5.5, (0.5, 1, 2, 3)), (-5.5, (0.5, 1, 2, 3))] ) -def test_generate_uv_offset_size(phase_centre, pixel_size): - phase_centre = phase_centre * apu.arcsec +def test_generate_uv_offset_size(phase_center, pixel_size): + phase_center = phase_center * apu.arcsec pixel_size = pixel_size * apu.arcsec / apu.pix m = 33 n = 32 @@ -93,18 +93,18 @@ def test_generate_uv_offset_size(phase_centre, pixel_size): for size in pixel_size: # Odd odd = np.linspace(-((m - 1) / 2) * (1 / (size * m)), ((m - 1) / 2) * (1 / (size * m)), m) / apu.pix - if phase_centre != 0: - odd += 1 / phase_centre + if phase_center != 0: + odd += 1 / phase_center assert np.allclose( - odd, generate_uv(m * apu.pix, phase_centre=phase_centre, pixel_size=size), atol=1e-8 / apu.arcsec + odd, generate_uv(m * apu.pix, phase_center=phase_center, pixel_size=size), atol=1e-8 / apu.arcsec ) # Even even = (np.arange(n) - n / 2 + 0.5) * (1 / (size * n)) / apu.pix - if phase_centre != 0: - even += 1 / phase_centre + if phase_center != 0: + even += 1 / phase_center assert np.allclose( - even, generate_uv(n * apu.pix, phase_centre=phase_centre, pixel_size=size), atol=1e-8 / apu.arcsec + even, generate_uv(n * apu.pix, phase_center=phase_center, pixel_size=size), atol=1e-8 / apu.arcsec ) @@ -114,10 +114,10 @@ def test_generate_uv_offset_size(phase_centre, pixel_size): def test_dft_idft(shape, pixel_size, center): data = np.arange(np.prod(shape)).reshape(shape) uu = generate_uv( - shape[0] * apu.pix, phase_centre=center[0] * apu.arcsec, pixel_size=pixel_size[0] * apu.arcsec / apu.pix + shape[0] * apu.pix, phase_center=center[0] * apu.arcsec, pixel_size=pixel_size[0] * apu.arcsec / apu.pix ) vv = generate_uv( - shape[1] * apu.pix, phase_centre=center[1] * apu.arcsec, pixel_size=pixel_size[1] * apu.arcsec / apu.pix + shape[1] * apu.pix, phase_center=center[1] * apu.arcsec, pixel_size=pixel_size[1] * apu.arcsec / apu.pix ) u, v = np.meshgrid(uu, vv, indexing="ij") u = u.flatten() @@ -213,45 +213,45 @@ def test_dft_idft(shape, pixel_size, center): # # # # @pytest.mark.skip('UV coordinate generation off somewhere') -# @pytest.mark.parametrize("phase_centre", [(0, 0), (2.1, 1.1), (5.4, -4.5), (-5.6, 5.6)]) -# def test_dft_idft_map_center(phase_centre): -# phase_centre = phase_centre * apu.arcsec +# @pytest.mark.parametrize("phase_center", [(0, 0), (2.1, 1.1), (5.4, -4.5), (-5.6, 5.6)]) +# def test_dft_idft_map_center(phase_center): +# phase_center = phase_center * apu.arcsec # shape = (33, 33) # m, n = shape # size = m * n # shape = shape * apu.pix -# uu = generate_uv(n * apu.pix, phase_centre=phase_centre[0]) -# vv = generate_uv(m * apu.pix, phase_centre=phase_centre[1]) +# uu = generate_uv(n * apu.pix, phase_center=phase_center[0]) +# vv = generate_uv(m * apu.pix, phase_center=phase_center[1]) # u, v = np.meshgrid(uu, vv) # u = u.flatten() # v = v.flatten() # # # All zeros # zeros = np.zeros((m, n)) -# vis = dft_map(zeros, u=u, v=v, phase_centre=phase_centre) +# vis = dft_map(zeros, u=u, v=v, phase_center=phase_center) # # All visibilities should be zero # assert_array_equal(vis, np.zeros(size, dtype=complex)) -# out_map = idft_map(vis, u=u, v=v, shape=shape, phase_centre=phase_centre) +# out_map = idft_map(vis, u=u, v=v, shape=shape, phase_center=phase_center) # # Should get back the original map after dft(idft()) # assert_array_equal(out_map / np.prod((m, n)), zeros) # # # All ones # ones = np.ones((m, n)) -# vis = dft_map(ones, u=u, v=v, phase_centre=phase_centre) -# out_map = idft_map(vis, u=u, v=v, shape=shape, phase_centre=phase_centre) +# vis = dft_map(ones, u=u, v=v, phase_center=phase_center) +# out_map = idft_map(vis, u=u, v=v, shape=shape, phase_center=phase_center) # assert_allclose(out_map / np.prod((m, n)), ones) # # # Delta # delta = zeros[:, :] # delta[m // 2, n // 2] = 1.0 -# vis = dft_map(delta, u=u, v=v, phase_centre=phase_centre) -# out_map = idft_map(vis, u=u, v=v, shape=shape, phase_centre=phase_centre) +# vis = dft_map(delta, u=u, v=v, phase_center=phase_center) +# out_map = idft_map(vis, u=u, v=v, shape=shape, phase_center=phase_center) # assert_allclose(out_map / np.prod((m, n)), delta, atol=1e-14) # # # Gaussian - astropy has axis in reverse order compared to numpy # gaussian = Gaussian2DKernel(5, x_size=n, y_size=m).array -# vis = dft_map(gaussian, u=u, v=v, phase_centre=phase_centre) -# out_map = idft_map(vis, u=u, v=v, shape=shape, phase_centre=phase_centre) +# vis = dft_map(gaussian, u=u, v=v, phase_center=phase_center) +# out_map = idft_map(vis, u=u, v=v, shape=shape, phase_center=phase_center) # assert_allclose(out_map / np.prod((m, n)), gaussian) # # @@ -337,15 +337,15 @@ def test_equivalence_of_convolve(): assert np.allclose(bp2, conv) assert np.allclose(bp2, bp3) - # Due to the enlarged psf need to only use the centre portion + # Due to the enlarged psf need to only use the center portion assert np.allclose(psf1[33:66, 33:66], psf2) -def test_phase_centre_equivalence(): +def test_phase_center_equivalence(): # Calculate full u, v coverage so will be equivalent to a discrete Fourier transform (DFT) data = np.random.randn(8, 8) / apu.arcsec**2 - u = generate_uv(8 * apu.pix, phase_centre=0 * apu.arcsec, pixel_size=1 * apu.arcsec / apu.pix) - v = generate_uv(8 * apu.pix, phase_centre=0 * apu.arcsec, pixel_size=1 * apu.arcsec / apu.pix) + u = generate_uv(8 * apu.pix, phase_center=0 * apu.arcsec, pixel_size=1 * apu.arcsec / apu.pix) + v = generate_uv(8 * apu.pix, phase_center=0 * apu.arcsec, pixel_size=1 * apu.arcsec / apu.pix) u, v = np.meshgrid(u, v) u, v = np.array([u, v]).reshape(2, u.size) / apu.arcsec @@ -360,11 +360,11 @@ def test_phase_centre_equivalence(): phase = np.arctan2(np.imag(vis), np.real(vis)) amp = np.abs(vis) - # change vis to a phase centre of 5, 5 + # change vis to a phase center of 5, 5 phase_shift = 2 * np.pi * (5 * apu.arcsec * u + 5 * apu.arcsec * v) * apu.rad vis_shifted = (np.cos(phase + phase_shift) + np.sin(phase + phase_shift) * 1j) * amp - # make image with centre of 5, 5 with shifted vis + # make image with center of 5, 5 with shifted vis img2 = idft_map( vis_shifted, u=u, @@ -372,7 +372,7 @@ def test_phase_centre_equivalence(): weights=1 / u.size, pixel_size=[1, 1] * apu.arcsec / apu.pix, shape=data.shape * apu.pix, - phase_centre=[5, 5] * apu.arcsec, + phase_center=[5, 5] * apu.arcsec, ) assert np.allclose(data, img2) @@ -387,16 +387,16 @@ def test_fft_equivalence(): data = np.arange(np.prod(shape)).reshape(shape) uu = generate_uv( - shape[0] * apu.pix, phase_centre=center[0] * apu.arcsec, pixel_size=pixel[0] * apu.arcsec / apu.pix + shape[0] * apu.pix, phase_center=center[0] * apu.arcsec, pixel_size=pixel[0] * apu.arcsec / apu.pix ) vv = generate_uv( - shape[1] * apu.pix, phase_centre=center[1] * apu.arcsec, pixel_size=pixel[1] * apu.arcsec / apu.pix + shape[1] * apu.pix, phase_center=center[1] * apu.arcsec, pixel_size=pixel[1] * apu.arcsec / apu.pix ) u, v = np.meshgrid(uu, vv, indexing="ij") u = u.flatten() v = v.flatten() - vis = dft_map(data, u=u, v=v, pixel_size=pixel * apu.arcsec / apu.pix, phase_centre=center * apu.arcsec) + vis = dft_map(data, u=u, v=v, pixel_size=pixel * apu.arcsec / apu.pix, phase_center=center * apu.arcsec) ft = fft2(data) fts = fftshift(ft) diff --git a/xrayvision/tests/test_visibility.py b/xrayvision/tests/test_visibility.py index 668863b..56d4fe8 100644 --- a/xrayvision/tests/test_visibility.py +++ b/xrayvision/tests/test_visibility.py @@ -3,10 +3,22 @@ from astropy.coordinates import get_body from astropy.tests.helper import assert_quantity_allclose from astropy.time import Time +import numpy as np from numpy.testing import assert_array_equal import xrayvision.visibility as vm -from xrayvision.visibility import Visibility +from xrayvision.visibility import Visibility, Visibilities + + +def test_visibilities(): + visibilities = np.array([1]) * apu.ct + u = np.array([1]) / apu.deg + v = np.array([1]) / apu.deg + vis = Visibilities(visibilities, u, v) + assert vis.visibilities == 1 * apu.ct + assert vis.u == 1 / apu.deg + assert vis.v == 1 / apu.deg + assert_array_equal([0, 0] * apu.arcsec, vis.phase_center) def test_visibility(): diff --git a/xrayvision/transform.py b/xrayvision/transform.py index d1ed6b7..1ff4847 100644 --- a/xrayvision/transform.py +++ b/xrayvision/transform.py @@ -22,17 +22,17 @@ def generate_xy( number_pixels: Quantity[apu.pix], *, - phase_centre: Optional[Quantity[apu.arcsec]] = 0.0 * apu.arcsec, + phase_center: Optional[Quantity[apu.arcsec]] = 0.0 * apu.arcsec, pixel_size: Optional[Quantity[apu.arcsec / apu.pix]] = 1.0 * apu.arcsec / apu.pix, ) -> Quantity[apu.arcsec]: """ - Generate the x or y coordinates given the number of pixels, phase_centre and pixel size. + Generate the x or y coordinates given the number of pixels, phase_center and pixel size. Parameters ---------- number_pixels : `int` Number of pixels - phase_centre : `float`, optional + phase_center : `float`, optional Center coordinates pixel_size : `float`, optional Size of pixel in physical units (e.g. arcsecs, meters) @@ -56,13 +56,13 @@ def generate_xy( >>> generate_xy(9*apu.pix, pixel_size=2.5 * apu.arcsec/apu.pix) - >>> generate_xy(9*apu.pix, phase_centre=10 * apu.arcsec, pixel_size=2.5 * apu.arcsec/apu.pix) + >>> generate_xy(9*apu.pix, phase_center=10 * apu.arcsec, pixel_size=2.5 * apu.arcsec/apu.pix) """ x = ( np.arange(number_pixels.to_value(apu.pixel)) - (number_pixels.to_value(apu.pix) / 2) + 0.5 - ) * apu.pix * pixel_size + phase_centre + ) * apu.pix * pixel_size + phase_center return x @@ -70,17 +70,17 @@ def generate_xy( def generate_uv( number_pixels: Quantity[apu.pix], *, - phase_centre: Optional[Quantity[apu.arcsec]] = 0.0 * apu.arcsec, + phase_center: Optional[Quantity[apu.arcsec]] = 0.0 * apu.arcsec, pixel_size: Optional[Quantity[apu.arcsec / apu.pix]] = 1.0 * apu.arcsec / apu.pix, ) -> Quantity[1 / apu.arcsec]: """ - Generate the u or v coordinates given the number of pixels, phase_centre and pixel size. + Generate the u or v coordinates given the number of pixels, phase_center and pixel size. Parameters ---------- number_pixels : `int` Number of pixels - phase_centre : `float`, optional + phase_center : `float`, optional Center coordinates pixel_size : `float`, optional Size of pixel in physical units (e.g. arcsecs, meters) @@ -106,7 +106,7 @@ def generate_uv( - >>> generate_uv(9*apu.pix, phase_centre=10 * apu.arcsec, pixel_size=2.5 * apu.arcsec/apu.pix) + >>> generate_uv(9*apu.pix, phase_center=10 * apu.arcsec, pixel_size=2.5 * apu.arcsec/apu.pix) @@ -117,8 +117,8 @@ def generate_uv( pixel_size * number_pixels ) - if phase_centre.value != 0.0: # type: ignore - x += 1 / phase_centre # type: ignore + if phase_center.value != 0.0: # type: ignore + x += 1 / phase_center # type: ignore return x @@ -128,7 +128,7 @@ def dft_map( *, u: Quantity[1 / apu.arcsec], v: Quantity[1 / apu.arcsec], - phase_centre: Quantity[apu.arcsec] = (0.0, 0.0) * apu.arcsec, + phase_center: Quantity[apu.arcsec] = (0.0, 0.0) * apu.arcsec, pixel_size: Quantity[apu.arcsec / apu.pix] = (1.0, 1.0) * apu.arcsec / apu.pix, ) -> Union[Quantity, npt.NDArray]: r""" @@ -142,8 +142,8 @@ def dft_map( Array of 2xN u coordinates where the visibilities are evaluated. v : Array of 2xN v coordinates where the visibilities are evaluated. - phase_centre : - Coordinates of the phase_centre of the map e.g. ``(0,0)`` or ``[5.0, -2.0]``. + phase_center : + Coordinates of the phase_center of the map e.g. ``(0,0)`` or ``[5.0, -2.0]``. pixel_size : `float` (dx, dy), optional The pixel size need not be square e.g. ``(1, 3)``. @@ -154,8 +154,8 @@ def dft_map( """ m, n = input_array.shape * apu.pix - x = generate_xy(m, phase_centre=phase_centre[0], pixel_size=pixel_size[0]) # type: ignore - y = generate_xy(n, phase_centre=phase_centre[1], pixel_size=pixel_size[1]) # type: ignore + x = generate_xy(m, phase_center=phase_center[0], pixel_size=pixel_size[0]) # type: ignore + y = generate_xy(n, phase_center=phase_center[1], pixel_size=pixel_size[1]) # type: ignore x, y = np.meshgrid(x, y, indexing="ij") uv = np.vstack([u, v]) @@ -190,7 +190,7 @@ def idft_map( v: Quantity[1 / apu.arcsec], shape: Quantity[apu.pix], weights: Optional[npt.NDArray] = None, - phase_centre: Quantity[apu.arcsec] = (0.0, 0.0) * apu.arcsec, + phase_center: Quantity[apu.arcsec] = (0.0, 0.0) * apu.arcsec, pixel_size: Quantity[apu.arcsec / apu.pix] = (1.0, 1.0) * apu.arcsec / apu.pix, ) -> Union[Quantity, npt.NDArray]: r""" @@ -208,8 +208,8 @@ def idft_map( The shape of the output array to create. weights : Array of weights for visibilities. - phase_centre : - Coordinates of the phase_centre of the map e.g. ``(0,0)`` or ``[5.0, -2.0]``. + phase_center : + Coordinates of the phase_center of the map e.g. ``(0,0)`` or ``[5.0, -2.0]``. pixel_size : The pixel size this need not be square e.g. ``(1, 3)``. @@ -220,8 +220,8 @@ def idft_map( """ m, n = shape - x = generate_xy(m, phase_centre=phase_centre[0], pixel_size=pixel_size[0]) # type: ignore - y = generate_xy(n, phase_centre=phase_centre[1], pixel_size=pixel_size[1]) # type: ignore + x = generate_xy(m, phase_center=phase_center[0], pixel_size=pixel_size[0]) # type: ignore + y = generate_xy(n, phase_center=phase_center[1], pixel_size=pixel_size[1]) # type: ignore x, y = np.meshgrid(x, y, indexing="ij") diff --git a/xrayvision/visibility.py b/xrayvision/visibility.py index f7e4dfe..4cd4ce8 100644 --- a/xrayvision/visibility.py +++ b/xrayvision/visibility.py @@ -523,7 +523,7 @@ def __init__( phase_centre: Phase centre of the visibility, defaults to (0,0). offset: - Offset of the phase_centre visibility, defaults to (0,0). + Offset of the phase_center visibility, defaults to (0,0). """ From 2f9bde916553a80dacb87551313dafbf6eddd0b6 Mon Sep 17 00:00:00 2001 From: paolomassa Date: Wed, 29 May 2024 12:58:08 +0200 Subject: [PATCH 04/12] Add changelog --- changelog/67.feature.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 changelog/67.feature.rst diff --git a/changelog/67.feature.rst b/changelog/67.feature.rst new file mode 100644 index 0000000..7326410 --- /dev/null +++ b/changelog/67.feature.rst @@ -0,0 +1 @@ +Make the software compatible with :class:`~xrayvision.visibility.Visibilities`. \ No newline at end of file From 3b30e1ee3498e5a9412ac6bc054a080f1b2e62f4 Mon Sep 17 00:00:00 2001 From: paolomassa Date: Wed, 29 May 2024 13:15:08 +0200 Subject: [PATCH 05/12] e --- changelog/67.feature.rst | 2 +- xrayvision/imaging.py | 21 +++++++++++++++++---- xrayvision/mem.py | 2 +- xrayvision/tests/test_imaging.py | 2 +- xrayvision/tests/test_visibility.py | 4 ++-- 5 files changed, 22 insertions(+), 9 deletions(-) diff --git a/changelog/67.feature.rst b/changelog/67.feature.rst index 7326410..b70d7af 100644 --- a/changelog/67.feature.rst +++ b/changelog/67.feature.rst @@ -1 +1 @@ -Make the software compatible with :class:`~xrayvision.visibility.Visibilities`. \ No newline at end of file +Make the software compatible with :class:`~xrayvision.visibility.Visibilities`. diff --git a/xrayvision/imaging.py b/xrayvision/imaging.py index d23de65..91f1a11 100644 --- a/xrayvision/imaging.py +++ b/xrayvision/imaging.py @@ -6,7 +6,7 @@ from sunpy.map import GenericMap, Map from xrayvision.transform import dft_map, idft_map -from xrayvision.visibility import Visibility, Visibilities +from xrayvision.visibility import Visibilities __all__ = [ "get_weights", @@ -161,7 +161,13 @@ def vis_to_image( shape = shape.to(apu.pixel) weights = get_weights(vis, scheme=scheme) bp_arr = idft_map( - vis.visibilities, u=vis.u, v=vis.v, shape=shape, weights=weights, pixel_size=pixel_size, phase_center=vis.phase_center + vis.visibilities, + u=vis.u, + v=vis.v, + shape=shape, + weights=weights, + pixel_size=pixel_size, + phase_center=vis.phase_center, ) return bp_arr @@ -237,7 +243,12 @@ def vis_psf_image( # Make sure psf is always odd so power is in exactly one pixel shape = [s // 2 * 2 + 1 for s in shape.to_value(apu.pix)] * shape.unit psf_arr = idft_map( - np.ones(vis.visibilities.shape) * vis.visibilities.unit, u=vis.u, v=vis.v, shape=shape, weights=weights, pixel_size=pixel_size + np.ones(vis.visibilities.shape) * vis.visibilities.unit, + u=vis.u, + v=vis.v, + shape=shape, + weights=weights, + pixel_size=pixel_size, ) return psf_arr @@ -346,6 +357,8 @@ def map_to_vis(amap: GenericMap, *, u: Quantity[1 / apu.arcsec], v: Quantity[1 / if "cdelt2" in meta: new_psize[0] = float(meta["cdelt2"]) - vis = image_to_vis(amap.quantity, u=u, v=v, pixel_size=new_psize * apu.arcsec / apu.pix, phase_center=new_pos * apu.arcsec) + vis = image_to_vis( + amap.quantity, u=u, v=v, pixel_size=new_psize * apu.arcsec / apu.pix, phase_center=new_pos * apu.arcsec + ) return vis diff --git a/xrayvision/mem.py b/xrayvision/mem.py index 3c68470..3dbdd98 100644 --- a/xrayvision/mem.py +++ b/xrayvision/mem.py @@ -27,7 +27,7 @@ "mem", ] -from xrayvision.visibility import Visibility, Visibilities +from xrayvision.visibility import Visibilities logger = get_logger(__name__, "DEBUG") diff --git a/xrayvision/tests/test_imaging.py b/xrayvision/tests/test_imaging.py index d19981f..0ac2eee 100644 --- a/xrayvision/tests/test_imaging.py +++ b/xrayvision/tests/test_imaging.py @@ -7,7 +7,7 @@ from xrayvision.imaging import image_to_vis, map_to_vis, vis_psf_image, vis_to_image, vis_to_map from xrayvision.transform import dft_map, generate_uv, idft_map -from xrayvision.visibility import Visibility, Visibilities +from xrayvision.visibility import Visibilities @pytest.fixture diff --git a/xrayvision/tests/test_visibility.py b/xrayvision/tests/test_visibility.py index 56d4fe8..037388b 100644 --- a/xrayvision/tests/test_visibility.py +++ b/xrayvision/tests/test_visibility.py @@ -1,13 +1,13 @@ import astropy.units as apu +import numpy as np import pytest from astropy.coordinates import get_body from astropy.tests.helper import assert_quantity_allclose from astropy.time import Time -import numpy as np from numpy.testing import assert_array_equal import xrayvision.visibility as vm -from xrayvision.visibility import Visibility, Visibilities +from xrayvision.visibility import Visibilities, Visibility def test_visibilities(): From 5d96f4142ebcca7b64f8da54a61d9239a7664a80 Mon Sep 17 00:00:00 2001 From: paolomassa Date: Fri, 31 May 2024 17:26:37 +0200 Subject: [PATCH 06/12] Fix STIX and RHESSI examples. Fix vis_to_image and image_to_vis phase center --- examples/rhessi.py | 40 ++++++++++--- examples/stix.py | 20 +++++-- xrayvision/imaging.py | 4 +- xrayvision/mem.py | 111 +++++++++++++++++++---------------- xrayvision/tests/test_mem.py | 2 +- 5 files changed, 108 insertions(+), 69 deletions(-) diff --git a/examples/rhessi.py b/examples/rhessi.py index bd75884..9906983 100644 --- a/examples/rhessi.py +++ b/examples/rhessi.py @@ -13,8 +13,8 @@ from xrayvision.clean import vis_clean from xrayvision.imaging import vis_psf_map, vis_to_map -from xrayvision.mem import mem -from xrayvision.visibility import Visibility +from xrayvision.mem import mem, resistant_mean +from xrayvision.visibility import Visibilities, VisMeta ############################################################################### # We will use `astropy.io.fits` to download and open the RHESSI visibility fits @@ -50,15 +50,19 @@ ############################################################################### # Now we can create the visibility object from the filtered visibilities. +meta = VisMeta({'vis_labels': vis_data["isc"]}) + vunit = apu.Unit("photon/(cm**2 s)") -vis = Visibility( - vis=vis_data["obsvis"] * vunit, - u=vis_data["u"] / apu.arcsec, - v=vis_data["v"] / apu.arcsec, - offset=vis_data["xyoffset"][0] * apu.arcsec, +vis = Visibilities( + vis_data["obsvis"] * vunit, + vis_data["u"] / apu.arcsec, + vis_data["v"] / apu.arcsec, + vis_data["xyoffset"][0] * apu.arcsec, + meta = meta, + amplitude_uncertainty = vis_data["sigamp"] * vunit ) -setattr(vis, "amplitude_error", vis_data["sigamp"] * vunit) -setattr(vis, "isc", vis_data["isc"]) +# setattr(vis, "amplitude_error", ) +# setattr(vis, "isc", vis_data["isc"]) ############################################################################### @@ -93,6 +97,24 @@ ############################################################################### # MEM +# Compute percent_lambda +# Loop through ISCs starting with 6-9, but if we don't have at least 2 vis, lower isc_min to include next one down, etc. +isc_min = 6 +nbig = 0 + +while isc_min >= 0 and nbig < 2: + ibig = np.argwhere(vis.meta.vis_labels >= isc_min) + nbig = len(ibig) + isc_min = isc_min - 1 + +# If still don't have at least 2 vis, return -1, otherwise calculate mean (but reject points > sigma away from mean) +if nbig < 2: + snr_value = -1 +else: + snr_value, _ = resistant_mean((np.abs(vis.visibilities[ibig]) / vis.amplitude_uncertainty[ibig]).flatten(), 3) + +percent_lambda = 11. / (snr_value**2 + 383.) + mem_map = mem(vis, shape=[129, 129] * apu.pixel, pixel_size=[2, 2] * apu.arcsec / apu.pix) mem_map.plot() diff --git a/examples/stix.py b/examples/stix.py index 5b908e6..8c54dc5 100644 --- a/examples/stix.py +++ b/examples/stix.py @@ -15,7 +15,9 @@ from xrayvision.clean import vis_clean from xrayvision.imaging import vis_psf_map, vis_to_map -from xrayvision.mem import mem +from xrayvision.mem import mem, resistant_mean + +import numpy as np ############################################################################### # Create images from STIX visibility data. @@ -23,11 +25,12 @@ # The STIX data has already been prepared and stored in python pickle format # the variables can be simply restored. -stix_data = pickle.load(urllib.request.urlopen("https://pub099.cs.technik.fhnw.ch/demo/stix_vis.pkl")) +with open('./stix_vis.pkl', 'rb') as file: + stix_data = pickle.load(file) -time_range, energy_range, offset, stix_vis = stix_data -stix_vis.phase_center = [0, 0] * apu.arcsec -stix_vis.offset = offset +time_range = stix_data['time_range'] +energy_range = stix_data['energy_range'] +stix_vis = stix_data['stix_visibilities'] ############################################################################### # Lets have a look at the point spread function (PSF) or dirty beam @@ -56,7 +59,12 @@ ############################################################################### # MEM -mem_map = mem(stix_vis, shape=[129, 129] * apu.pixel, pixel_size=[2, 2] * apu.arcsec / apu.pix) +# Compute percent_lambda +snr_value, _ = resistant_mean((np.abs(stix_vis.visibilities) / stix_vis.amplitude_uncertainty).flatten(), 3) +percent_lambda = 2 / (snr_value**2 + 90) + +mem_map = mem(stix_vis, shape=[129, 129] * apu.pixel, pixel_size=[2, 2] * apu.arcsec / apu.pix, + percent_lambda=percent_lambda) mem_map.plot() ############################################################################### diff --git a/xrayvision/imaging.py b/xrayvision/imaging.py index 91f1a11..124fb32 100644 --- a/xrayvision/imaging.py +++ b/xrayvision/imaging.py @@ -125,7 +125,7 @@ def image_to_vis( pixel_size = validate_and_expand_kwarg(pixel_size, "pixel_size") if not (apu.get_physical_type((1 / u).unit) == ANGLE and apu.get_physical_type((1 / v).unit) == ANGLE): raise ValueError("u and v must be inverse angle (e.g. 1/deg or 1/arcsec") - vis = dft_map(image, u=u, v=v, phase_center=phase_center, pixel_size=pixel_size) + vis = dft_map(image, u=u, v=v, phase_center=[0., 0.]*apu.arcsec, pixel_size=pixel_size) # TODO: adapt to generic map center return Visibilities(vis, u=u, v=v, phase_center=phase_center) @@ -167,7 +167,7 @@ def vis_to_image( shape=shape, weights=weights, pixel_size=pixel_size, - phase_center=vis.phase_center, + phase_center=[0., 0.] * apu.arcsec, # TODO update to have generic image center ) return bp_arr diff --git a/xrayvision/mem.py b/xrayvision/mem.py index 3dbdd98..51cd1df 100644 --- a/xrayvision/mem.py +++ b/xrayvision/mem.py @@ -197,8 +197,12 @@ def _prepare_for_optimise(pixel, shape, vis): # 'Visib' is composed by the real and imaginary part of the visibilities Visib = np.concatenate([ReV, ImV], axis=-1) # Standard deviation of the real and imaginary part of the visibilities - sigma_Re = vis.amplitude_error - sigma_Im = vis.amplitude_error + if vis.amplitude_uncertainty is None: + sigma_Re = np.ones_like(ReV) + sigma_Im = np.ones_like(ImV) + else: + sigma_Re = vis.amplitude_uncertainty + sigma_Im = vis.amplitude_uncertainty # 'sigma': standard deviation of the data contained in 'Visib' sigma = np.concatenate([sigma_Re, sigma_Im], axis=-1) # RESCALING of 'Hv' AND 'Visib'(NEEDED FOR COMPUTING THE VALUE OF THE \chi ** 2; FUNCTION) @@ -234,6 +238,11 @@ def _get_mean_visibilities(vis, shape, pixel): Mean Visibilities """ + if vis.amplitude_uncertainty is None: + amplitude_uncertainty = np.ones_like(vis.visibilities) + else: + amplitude_uncertainty = vis.amplitude_uncertainty + imsize2 = shape[0] / 2 pixel_size = 1 / (shape[0] * pixel[0]) @@ -255,7 +264,7 @@ def _get_mean_visibilities(vis, shape, pixel): u = np.zeros(n_vis) * (1 / apu.arcsec) v = np.zeros(n_vis) * (1 / apu.arcsec) den = np.zeros(n_vis) - weights = np.zeros_like(vis.amplitude_error**2) + weights = np.ones_like(vis.visibilities**2) visib = np.zeros_like(vis.visibilities) for ip in range(n_vis): # what about 0.5 pix offset @@ -270,7 +279,7 @@ def _get_mean_visibilities(vis, shape, pixel): # to the position (i, j) of the discretization of the (u,v)-plane 'iuarr' visib[count] = vis.visibilities[ip] - weights[count] = vis.amplitude_error[ip] ** 2.0 + weights[count] = amplitude_uncertainty[ip] ** 2.0 den[count] = 1.0 iuarr[i, j] = count @@ -282,7 +291,7 @@ def _get_mean_visibilities(vis, shape, pixel): den[iuarr[i, j].astype(int)] += 1.0 # Save the sum of the variances of the amplitudes of the visibilities that # correspond to the same position (i, j) - weights[iuarr[i, j].astype(int)] += vis.amplitude_error[ip] ** 2 + weights[iuarr[i, j].astype(int)] += amplitude_uncertainty[ip] ** 2 u = u[:count] v = v[:count] @@ -296,7 +305,7 @@ def _get_mean_visibilities(vis, shape, pixel): # correspond to the same position in the discretization of the (u,v)-plane weights = np.sqrt(weights[:count]) / den - return SimpleNamespace(u=u, v=v, visibilities=visib, amplitude_error=weights) + return SimpleNamespace(u=u, v=v, visibilities=visib, amplitude_uncertainty=weights) def _proximal_entropy(y, m, lamba, Lip, tol=10**-10): @@ -565,48 +574,48 @@ def resistant_mean(data, sigma_cut): return mean, sigma -def _get_percent_lambda(vis): - r""" - Return 'percent_lambda' use with MEM - - Calculate the signal-to-noise ratio (SNR) for the given visibility bag by trying to use the - coarser sub-collimators adding finer ones until there are at least 2 visibilities, then use - resistant mean of of abs(obsvis) / sigamp - - Parameters - ---------- - vis - - Returns - ------- - - """ - # Loop through ISCs starting with 3-10, but if we don't have at least 2 vis, lower isc_min to - # include next one down, etc. - - # TODO this start at 3 not 10? - isc_min = 3 - nbig = 0 - - if hasattr(vis, "label"): - isc_sizes = np.array([float(s[:-1]) for s in vis.label]) - while isc_min >= 0 and nbig < 2: - ibig = np.argwhere(isc_sizes >= isc_min) - isc_min = isc_min - 1 - elif hasattr(vis, "isc"): - ibig = np.arange(vis.isc.size) - - # If still don't have at least 2 vis, return -1, otherwise calculate mean - # (but reject points > sigma away from mean) - if ibig.size < 2: - snr_value = -1 - else: - snr_value, _ = resistant_mean((np.abs(vis.visibilities[ibig]) / vis.amplitude_error[ibig]).flatten(), 3) - - # TODO magic numbers - percent_lambda = 2 / (snr_value**2 + 90) - - return percent_lambda +# def _get_percent_lambda(vis): +# r""" +# Return 'percent_lambda' use with MEM +# +# Calculate the signal-to-noise ratio (SNR) for the given visibility bag by trying to use the +# coarser sub-collimators adding finer ones until there are at least 2 visibilities, then use +# resistant mean of of abs(obsvis) / sigamp +# +# Parameters +# ---------- +# vis +# +# Returns +# ------- +# +# """ +# # Loop through ISCs starting with 3-10, but if we don't have at least 2 vis, lower isc_min to +# # include next one down, etc. +# +# # TODO this start at 3 not 10? +# isc_min = 3 +# nbig = 0 +# +# if hasattr(vis.meta, "label"): +# isc_sizes = np.array([float(s[:-1]) for s in vis.meta.label]) +# while isc_min >= 0 and nbig < 2: +# ibig = np.argwhere(isc_sizes >= isc_min) +# isc_min = isc_min - 1 +# elif hasattr(vis.meta, "isc"): +# ibig = np.arange(vis.meta.isc.size) +# +# # If still don't have at least 2 vis, return -1, otherwise calculate mean +# # (but reject points > sigma away from mean) +# if ibig.size < 2: +# snr_value = -1 +# else: +# snr_value, _ = resistant_mean((np.abs(vis.visibilities[ibig]) / vis.amplitude_error[ibig]).flatten(), 3) +# +# # TODO magic numbers +# percent_lambda = 2 / (snr_value**2 + 90) +# +# return percent_lambda @apu.quantity_input @@ -615,7 +624,7 @@ def mem( shape: Quantity[apu.pix], pixel_size: Quantity[apu.arcsec / apu.pix], *, - percent_lambda: Optional[Quantity[apu.percent]] = None, + percent_lambda: Optional[Quantity[apu.percent]] = 0.02 * apu.percent, maxiter: int = 1000, tol: float = 1e-3, map: bool = True, @@ -645,8 +654,8 @@ def mem( """ total_flux = _estimate_flux(vis, shape, pixel_size) - if percent_lambda is None: - percent_lambda = _get_percent_lambda(vis) + # if percent_lambda is None: + # percent_lambda = _get_percent_lambda(vis) mean_vis = _get_mean_visibilities(vis, shape, pixel_size) Hv, Lip, Visib = _prepare_for_optimise(pixel_size, shape, mean_vis) diff --git a/xrayvision/tests/test_mem.py b/xrayvision/tests/test_mem.py index 69a41ee..d1079d1 100644 --- a/xrayvision/tests/test_mem.py +++ b/xrayvision/tests/test_mem.py @@ -61,6 +61,6 @@ def test_mem(): setattr(vis, "label", [str(x) for x in np.sqrt(x**2 + y**2).flatten()]) res = mem( - vis, percent_lambda=None, shape=(m, n) * u.pix, pixel_size=[2, 2] * u.arcsec / u.pix, maxiter=1000, tol=1e-3 + vis, shape=(m, n) * u.pix, pixel_size=[2, 2] * u.arcsec / u.pix, maxiter=1000, tol=1e-3 ) assert_allclose(res.data, data, atol=1e-1) From 42e85fc093a240fb08e247e0a1cd73dc05b6559c Mon Sep 17 00:00:00 2001 From: paolomassa Date: Fri, 31 May 2024 17:29:46 +0200 Subject: [PATCH 07/12] Add STIX visibility file. --- examples/stix_vis.pkl | Bin 0 -> 4633 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 examples/stix_vis.pkl diff --git a/examples/stix_vis.pkl b/examples/stix_vis.pkl new file mode 100644 index 0000000000000000000000000000000000000000..d44eeb77c1c5f29e35b47c647028c7914b9f8e61 GIT binary patch literal 4633 zcmbtY3s4kC8is{k-lzzO7<2_6pu9n%XzW#tLfE&1@cedKHsgnM-`s9hzti!3P>OCMos)p4sIk$z9z+{nI@? zfB*gckN*4bUG2LmG|Wr=tC!=YVNQXqrIxc!ww$%H0*jrp6j965b!5oSCHH)(lcnI2 zzHF^>!tZ1v!k^LMhtjfR#iG);=qVCIT+gM7Nl4EvNX&+*jKmgd2N$!0r`lgW@13m{>RyQdnikXX(@ z(SpD`i4~bg$irU>qn!e8chn`TGF1_>6p~}NIqXhWz>8p2xySMbsvKlHfWbEM60mDI*)y`G+j4~NC*s8+L+h{?R<68lCl9>$2*IvmwNXM`omI5`H3=vo< zi``67U{Mu-G2SdB$YMF5_3f`~5NPM%u@iGd+yp6bK94Kd+I~#T6HL4ol+5H(nMr0+ ziGJQTP!GI_`Cj#yM~Tgfm}1OK1I{9vnGQ`c%{*2wGqWY5%ARc`<~0(qMtgE!_9*i_^E+WbmwKJ#ny^V~-2Z?4Dx%hA$Y)U72D zV&!<|;dui&+>X#ySIUxqHEXq3^0aNBGa8si%^3}v28|vKVm&hk_`^m!otih2NyIGd z@i(L}gA)wTSq84Lf!V07u|Zp7qh}4H#y?;M;+ftRF!Kzko9~_Pop=Aig9n&d%B=D@ z)4LGXEW|Z}or4aib+o0ns!PQ)FEHD%JG??P&TJ#>6Ejx9_X;O+JyAT+kPDe6rp3*w z>w{yLSl)vq528KKEQj5Vsq4dJ=sU|pvl0X^+i=tt1d9XwdynuS;`w6`%V$U`(!tdjR|2&AT2YL5MgsBQd|3WkjvhPUJdw4hY zNrgTsA-dqm{EhKES-2(Bi``~q1gmj*+$rxOL zw3LjjRH8Jv_(F3IPs)o-f!$AEyI3^j-oJYdz(q5FZH z6_9kynfB_+J|HzNf>UIo7@gz?WnBR&a@4vaUJ4-PoRu!eRw&3g(FVncTFKZ%tH7$2 z4Qy#rYuGucRYa&BT4gERQ@02-6df?ZD%waY%+Wfna=urSoa0?Cllp=r(WoZKPpWii zEkd0G6EQkAQ(mvuAo2M=lJ^X`u&_|u5M)>6dkJ#+`kXQ!NC)9Skw8O%LV=7x5kS#E zzCckx{y>ia4FMVjGzQ23G#W?;G!jS;Gy=#ECfR#iNq7*2ly{>!C6R<#G&6;KJ2x>^=Ru$0+u=Pm96sTt18a3->a)QYT) zRDMi!ev}g`PRv2fvRtKH2!({%YIkDj`BN-s!L?9^L0kt?%W0k^=G7`je<{3gq3QVu zM{J3kPxB}#HLJ#hUQBe8%*0grLAPosJS0EPqWM#EPgO7s_VG_u7-8};r)Bl}q$(1| z264*l)iB<9=2Y3NXW~s|)iGyc%imWf&i9p$Zc8;yplY|z4>hPSHY;lnb#(nT3$~pySwY3=+uEr zt|`CmY$}>LFTu2WZ06yU#yd({>xihV@&wcA(J$-y*;kOGW5<-a{t2et$8IF+fAV#O zxfYc;-XCvrjTE=PRp#ru(>6<&sUL5e{~+Y=b+fK0#ZfzsG8N-Zr)=k?%RV~ODtcAQ zwfnKA#F!H$-M9TrU!0*nyv&a_g*P`ozv8i5Xy3slKd`Sxn-Wjm*tF)9!Iif7*q2*V zB1|*>w7K1$a9_Fd>`x(MA0KMkz5Dsg?b$k0>9uW^r_+2*$G1-TzAntqCI04kVquEj zl=zn~K8?L{3Qd}}^G@yl>u7l3-o&kKHpawPhduany(#=+4=;w#K4vl(}nq zh_XLFWX0xWbfM{oX2&!qb)gRyH$46EU%x|M<2gf2S|@sE#oYQpughrIxz~@rAzenF zzWx1y;X5y(_|}roLKc66&TW0iaR25-)O`2Y2MyM*(AjV_30dk^Tsm^`JeqoG z0##6S9*vV6hgvUO7{qAd|C`a~fi>6XxDu~(oyR{hqsZ%b*Bx0HFeAzRKBaZh*sQS* z)YY+O)S>DU<=(iCEPg0fOdlD(9X+B38O&dQcVRq-TKG@i(+8EIm&3a2qleU^m)tE| zYJKipW9;EufBrZ(@Z2h;v*e)v!QE?>pVvP-E~?t0gck0`F{^6O*~q2^k)16Hz3xzK z$?!F3$BMh(-inw$qdll#Vq4B8^qRZ9efX_OI|A0Bj+zbCjhi4TT z{J-q6Y|V$%lI{{^5pGC1eWwCBw0G1IKkLho6LXc1*z+@|p5{>f(#xhd<5*?#4;wbb zJvkrk?K&}O@=wc@*XFUq?}RKzhIVm}+@g--q(gW_6H-Hl~*!MREWwa8buJoN{q=D~^sksG(5w&EjG4I9^? zom#6*Y>wJ{x?!iX@#!SJ(Di~c>9PHJ>zUez<3l7nOA2xS%0RAXaI;jisubxt1<9S! znY>J-;dV!U_3$t9bO%YbI@p0)@Q}mI%geIWrMBmrT;L4I{Ce#Y&zy9^$|9UJ3VSrDi3iKDsbHEjwdcc513Ta5OVLJ9C;0_js+^@134+$PdI~^1lJV C&>ph@ literal 0 HcmV?d00001 From e857a0666e67db2c5c176406bc0787a029f377f9 Mon Sep 17 00:00:00 2001 From: paolomassa Date: Fri, 31 May 2024 17:39:45 +0200 Subject: [PATCH 08/12] Cleaup mem.py --- xrayvision/mem.py | 44 -------------------------------------------- 1 file changed, 44 deletions(-) diff --git a/xrayvision/mem.py b/xrayvision/mem.py index 51cd1df..422d5f5 100644 --- a/xrayvision/mem.py +++ b/xrayvision/mem.py @@ -574,50 +574,6 @@ def resistant_mean(data, sigma_cut): return mean, sigma -# def _get_percent_lambda(vis): -# r""" -# Return 'percent_lambda' use with MEM -# -# Calculate the signal-to-noise ratio (SNR) for the given visibility bag by trying to use the -# coarser sub-collimators adding finer ones until there are at least 2 visibilities, then use -# resistant mean of of abs(obsvis) / sigamp -# -# Parameters -# ---------- -# vis -# -# Returns -# ------- -# -# """ -# # Loop through ISCs starting with 3-10, but if we don't have at least 2 vis, lower isc_min to -# # include next one down, etc. -# -# # TODO this start at 3 not 10? -# isc_min = 3 -# nbig = 0 -# -# if hasattr(vis.meta, "label"): -# isc_sizes = np.array([float(s[:-1]) for s in vis.meta.label]) -# while isc_min >= 0 and nbig < 2: -# ibig = np.argwhere(isc_sizes >= isc_min) -# isc_min = isc_min - 1 -# elif hasattr(vis.meta, "isc"): -# ibig = np.arange(vis.meta.isc.size) -# -# # If still don't have at least 2 vis, return -1, otherwise calculate mean -# # (but reject points > sigma away from mean) -# if ibig.size < 2: -# snr_value = -1 -# else: -# snr_value, _ = resistant_mean((np.abs(vis.visibilities[ibig]) / vis.amplitude_error[ibig]).flatten(), 3) -# -# # TODO magic numbers -# percent_lambda = 2 / (snr_value**2 + 90) -# -# return percent_lambda - - @apu.quantity_input def mem( vis: Visibilities, From 261c7143503bbddc018f8f601156623c0e4e0e66 Mon Sep 17 00:00:00 2001 From: paolomassa Date: Fri, 31 May 2024 18:26:10 +0200 Subject: [PATCH 09/12] Fix pre-commit. --- examples/rhessi.py | 8 ++++---- examples/stix.py | 17 ++++++++--------- xrayvision/imaging.py | 6 ++++-- xrayvision/tests/test_mem.py | 4 +--- 4 files changed, 17 insertions(+), 18 deletions(-) diff --git a/examples/rhessi.py b/examples/rhessi.py index 9906983..be1dc8f 100644 --- a/examples/rhessi.py +++ b/examples/rhessi.py @@ -50,7 +50,7 @@ ############################################################################### # Now we can create the visibility object from the filtered visibilities. -meta = VisMeta({'vis_labels': vis_data["isc"]}) +meta = VisMeta({"vis_labels": vis_data["isc"]}) vunit = apu.Unit("photon/(cm**2 s)") vis = Visibilities( @@ -58,8 +58,8 @@ vis_data["u"] / apu.arcsec, vis_data["v"] / apu.arcsec, vis_data["xyoffset"][0] * apu.arcsec, - meta = meta, - amplitude_uncertainty = vis_data["sigamp"] * vunit + meta=meta, + amplitude_uncertainty=vis_data["sigamp"] * vunit, ) # setattr(vis, "amplitude_error", ) # setattr(vis, "isc", vis_data["isc"]) @@ -113,7 +113,7 @@ else: snr_value, _ = resistant_mean((np.abs(vis.visibilities[ibig]) / vis.amplitude_uncertainty[ibig]).flatten(), 3) -percent_lambda = 11. / (snr_value**2 + 383.) +percent_lambda = 11.0 / (snr_value**2 + 383.0) mem_map = mem(vis, shape=[129, 129] * apu.pixel, pixel_size=[2, 2] * apu.arcsec / apu.pix) mem_map.plot() diff --git a/examples/stix.py b/examples/stix.py index 8c54dc5..e4a36a3 100644 --- a/examples/stix.py +++ b/examples/stix.py @@ -8,29 +8,27 @@ """ import pickle -import urllib.request import astropy.units as apu import matplotlib.pyplot as plt +import numpy as np from xrayvision.clean import vis_clean from xrayvision.imaging import vis_psf_map, vis_to_map from xrayvision.mem import mem, resistant_mean -import numpy as np - ############################################################################### # Create images from STIX visibility data. # # The STIX data has already been prepared and stored in python pickle format # the variables can be simply restored. -with open('./stix_vis.pkl', 'rb') as file: +with open("./stix_vis.pkl", "rb") as file: stix_data = pickle.load(file) -time_range = stix_data['time_range'] -energy_range = stix_data['energy_range'] -stix_vis = stix_data['stix_visibilities'] +time_range = stix_data["time_range"] +energy_range = stix_data["energy_range"] +stix_vis = stix_data["stix_visibilities"] ############################################################################### # Lets have a look at the point spread function (PSF) or dirty beam @@ -63,8 +61,9 @@ snr_value, _ = resistant_mean((np.abs(stix_vis.visibilities) / stix_vis.amplitude_uncertainty).flatten(), 3) percent_lambda = 2 / (snr_value**2 + 90) -mem_map = mem(stix_vis, shape=[129, 129] * apu.pixel, pixel_size=[2, 2] * apu.arcsec / apu.pix, - percent_lambda=percent_lambda) +mem_map = mem( + stix_vis, shape=[129, 129] * apu.pixel, pixel_size=[2, 2] * apu.arcsec / apu.pix, percent_lambda=percent_lambda +) mem_map.plot() ############################################################################### diff --git a/xrayvision/imaging.py b/xrayvision/imaging.py index 124fb32..a185d52 100644 --- a/xrayvision/imaging.py +++ b/xrayvision/imaging.py @@ -125,7 +125,9 @@ def image_to_vis( pixel_size = validate_and_expand_kwarg(pixel_size, "pixel_size") if not (apu.get_physical_type((1 / u).unit) == ANGLE and apu.get_physical_type((1 / v).unit) == ANGLE): raise ValueError("u and v must be inverse angle (e.g. 1/deg or 1/arcsec") - vis = dft_map(image, u=u, v=v, phase_center=[0., 0.]*apu.arcsec, pixel_size=pixel_size) # TODO: adapt to generic map center + vis = dft_map( + image, u=u, v=v, phase_center=[0.0, 0.0] * apu.arcsec, pixel_size=pixel_size + ) # TODO: adapt to generic map center return Visibilities(vis, u=u, v=v, phase_center=phase_center) @@ -167,7 +169,7 @@ def vis_to_image( shape=shape, weights=weights, pixel_size=pixel_size, - phase_center=[0., 0.] * apu.arcsec, # TODO update to have generic image center + phase_center=[0.0, 0.0] * apu.arcsec, # TODO update to have generic image center ) return bp_arr diff --git a/xrayvision/tests/test_mem.py b/xrayvision/tests/test_mem.py index d1079d1..9c24ad6 100644 --- a/xrayvision/tests/test_mem.py +++ b/xrayvision/tests/test_mem.py @@ -60,7 +60,5 @@ def test_mem(): setattr(vis, "amplitude_error", np.sqrt(np.abs(vis.visibilities))) setattr(vis, "label", [str(x) for x in np.sqrt(x**2 + y**2).flatten()]) - res = mem( - vis, shape=(m, n) * u.pix, pixel_size=[2, 2] * u.arcsec / u.pix, maxiter=1000, tol=1e-3 - ) + res = mem(vis, shape=(m, n) * u.pix, pixel_size=[2, 2] * u.arcsec / u.pix, maxiter=1000, tol=1e-3) assert_allclose(res.data, data, atol=1e-1) From 3238c149efc992034c89b8ebef4035e858856c20 Mon Sep 17 00:00:00 2001 From: Paolo Massa <64010598+paolomassa@users.noreply.github.com> Date: Fri, 31 May 2024 21:39:42 +0200 Subject: [PATCH 10/12] Update examples/rhessi.py Co-authored-by: Shane Maloney --- examples/rhessi.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/examples/rhessi.py b/examples/rhessi.py index be1dc8f..a170dbb 100644 --- a/examples/rhessi.py +++ b/examples/rhessi.py @@ -61,8 +61,6 @@ meta=meta, amplitude_uncertainty=vis_data["sigamp"] * vunit, ) -# setattr(vis, "amplitude_error", ) -# setattr(vis, "isc", vis_data["isc"]) ############################################################################### From d55568f0eabbe36eaab3c9e787e0a4ea41533849 Mon Sep 17 00:00:00 2001 From: Paolo Massa <64010598+paolomassa@users.noreply.github.com> Date: Fri, 31 May 2024 21:40:26 +0200 Subject: [PATCH 11/12] Update xrayvision/mem.py Co-authored-by: Shane Maloney --- xrayvision/mem.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/xrayvision/mem.py b/xrayvision/mem.py index 422d5f5..0804213 100644 --- a/xrayvision/mem.py +++ b/xrayvision/mem.py @@ -610,8 +610,6 @@ def mem( """ total_flux = _estimate_flux(vis, shape, pixel_size) - # if percent_lambda is None: - # percent_lambda = _get_percent_lambda(vis) mean_vis = _get_mean_visibilities(vis, shape, pixel_size) Hv, Lip, Visib = _prepare_for_optimise(pixel_size, shape, mean_vis) From 900e12e4feffe9c2e497343bd4f8501298bdc55c Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Sat, 1 Jun 2024 12:20:45 +0100 Subject: [PATCH 12/12] Tiny updates --- examples/rhessi.py | 8 ++++---- xrayvision/tests/test_imaging.py | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/rhessi.py b/examples/rhessi.py index a170dbb..031dbc3 100644 --- a/examples/rhessi.py +++ b/examples/rhessi.py @@ -54,10 +54,10 @@ vunit = apu.Unit("photon/(cm**2 s)") vis = Visibilities( - vis_data["obsvis"] * vunit, - vis_data["u"] / apu.arcsec, - vis_data["v"] / apu.arcsec, - vis_data["xyoffset"][0] * apu.arcsec, + visibilities=vis_data["obsvis"] * vunit, + u=vis_data["u"] / apu.arcsec, + v=vis_data["v"] / apu.arcsec, + phase_center=vis_data["xyoffset"][0] * apu.arcsec, meta=meta, amplitude_uncertainty=vis_data["sigamp"] * vunit, ) diff --git a/xrayvision/tests/test_imaging.py b/xrayvision/tests/test_imaging.py index 0ac2eee..a359954 100644 --- a/xrayvision/tests/test_imaging.py +++ b/xrayvision/tests/test_imaging.py @@ -93,7 +93,7 @@ def test_vis_to_psf(pixel_size, uv): weights = np.sqrt(u**2 + v**2).value weights /= weights.sum() psf_calc = idft_map(obs_vis, u=u, v=v, shape=[65, 65] * apu.pix, pixel_size=ps, weights=weights) - vis = Visibilities(obs_vis, u, v) + vis = Visibilities(obs_vis, u=u, v=v) res = vis_psf_image(vis, shape=[65, 65] * apu.pixel, pixel_size=ps, scheme="uniform") assert_allclose(res, psf_calc)