Noise (& beam) response¶
Download¶
You can download the data directory data_tutorials using:
wget https://www.mso.anu.edu.au/~amarchal/shared/ivis/data_tutorials.zip
unzip data_tutorials.zip -d data_tutorials
rm data_tutorials.zip
Kindly provided by Enrico Di Teodoro (UniFI) and Karlie Noon (RSAA/ANU).
Method¶
To characterize the image-domain response of the noise in the visibility data, we performed imaging on a single realization constructed using the
SIGMAarray available in the Measurement Set.
We then computed the one-dimensional spatial power spectrum P(k) of the resulting image.Note that to fully characterize P(k), we recommend repeating this procedure over multiple random realizations (typically at least 20) and averaging the resulting spectra.
Read data¶
import glob
from tqdm import tqdm as tqdm
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import ScalarFormatter, FormatStrFormatter
from astropy.io import fits
from radio_beam import Beam
from astropy import units as u
from reproject import reproject_interp
from ivis.io import DataProcessor
from ivis.imager import Imager3D
from ivis.logger import logger
from ivis.utils import dutils, fourier
from ivis.readers import CasacoreReader
from ivis.types import VisIData
from ivis.models import ClassicIViS3D
path_ms = "./data_tutorials/ivis_data/msl_mw/" #directory of measurement sets
path_beams = "./data_tutorials/ivis_data/BEAMS/" #directory of primary beams
path_sd = None #path single-dish data
pathout = "./data_tutorials/ivis_data/" #path where data will be packaged and stored
#REF WCS INPUT USER
filename = "./data_tutorials/ivis_data/MW-C10_mom0th_NHI.fits"
target_header = fits.open(filename)[0].header
shape = (target_header["NAXIS2"],target_header["NAXIS1"])
#create data processor
data_processor = DataProcessor(path_ms, path_beams, path_sd, pathout)
# pre-compute pb and interpolation grids — this can be commented after first compute
logger.disabled = True
data_processor.compute_pb_and_grid(target_header, fitsname_pb="reproj_pb.fits", fitsname_grid="grid_interp.fits")
logger.disabled = False
pb, grid = data_processor.read_pb_and_grid(fitsname_pb="reproj_pb.fits", fitsname_grid="grid_interp.fits")
#Dummy sd array
sd = np.zeros(shape)
#Dummy Beam sd
beam_sd = Beam(1*u.deg, 1*u.deg, 1.e-12*u.deg)
[2025-08-25 07:25:51 UTC] [IViS] [info] [Initialize DataProcessor ]
[2025-08-25 07:25:51 UTC] [IViS] [info] -------------------------------------------------------------------------
[2025-08-25 07:25:51 UTC] [IViS] [info] Timestamp: 2025-08-25 07:25:51 UTC
[2025-08-25 07:25:51 UTC] [IViS] [info]
[2025-08-25 07:25:51 UTC] [IViS] [info] _ _| \ \ / _) ___|
[2025-08-25 07:25:51 UTC] [IViS] [info] | \ \ / | \___ \
[2025-08-25 07:25:51 UTC] [IViS] [info] | \ \ / | |
[2025-08-25 07:25:51 UTC] [IViS] [info] ___| \_/ _| _____/
[2025-08-25 07:25:51 UTC] [IViS] [info]
[2025-08-25 07:25:51 UTC] [IViS] [info] Version 1.0.0
[2025-08-25 07:25:51 UTC] [IViS] [info] IViS is released as open-source software
[2025-08-25 07:25:51 UTC] [IViS] [info] Author: @amarchal
[2025-08-25 07:25:51 UTC] [IViS] [info] Documentation: https://ivis-dev.readthedocs.io/en/latest/
[2025-08-25 07:25:51 UTC] [IViS] [info] Github: https://github.com/antoinemarchal/ivis
[2025-08-25 07:25:51 UTC] [IViS] [info] -------------------------------------------------------------------------
100%|█████████████████████████████████████████████| 5/5 [00:02<00:00, 1.96it/s]
# -------------------
# Read visibilities into VisIData dataclass
# -------------------
reader = CasacoreReader(
prefer_weight_spectrum=False,
keep_autocorr=False,
n_workers=4)
I: VisIData = reader.read_blocks_I(
path_ms,
uvmin=0.0,
uvmax=np.inf,
chan_sel=slice(0, 1),
)
[2025-08-25 07:25:55 UTC] [IViS] [info] [BLOCK] Loading block from: ./data_tutorials/ivis_data/msl_mw/
[2025-08-25 07:25:55 UTC] [IViS] [info] [BLOCK] Loading 5 beam(s) from: ./data_tutorials/ivis_data/msl_mw/
[2025-08-25 07:25:55 UTC] [IViS] [info] [BLOCK] Parallel read with 4 workers (order-preserving; selective)
[2025-08-25 07:25:55 UTC] [IViS] [info] [MS] Opening: ./data_tutorials/ivis_data/msl_mw/MW-C10_1_MW_chan_-32kms.ms
[2025-08-25 07:25:55 UTC] [IViS] [info] [MS] Opening: ./data_tutorials/ivis_data/msl_mw/MW-C10_3_MW_chan_-32kms.ms
[2025-08-25 07:25:55 UTC] [IViS] [info] [MS] Opening: ./data_tutorials/ivis_data/msl_mw/MW-C10_2_MW_chan_-32kms.ms
[2025-08-25 07:25:55 UTC] [IViS] [info] [MS] Opening: ./data_tutorials/ivis_data/msl_mw/MW-C10_4_MW_chan_-32kms.ms
[2025-08-25 07:25:56 UTC] [IViS] [info] [MS] Done: MW-C10_1_MW_chan_-32kms.ms rows=626313
[2025-08-25 07:25:56 UTC] [IViS] [info] [MS] Opening: ./data_tutorials/ivis_data/msl_mw/MW-C10_5_MW_chan_-32kms.ms
[2025-08-25 07:25:56 UTC] [IViS] [info] [MS] Done: MW-C10_4_MW_chan_-32kms.ms rows=1140552
[2025-08-25 07:25:56 UTC] [IViS] [info] [MS] Done: MW-C10_3_MW_chan_-32kms.ms rows=1140552
[2025-08-25 07:25:56 UTC] [IViS] [info] [MS] Done: MW-C10_2_MW_chan_-32kms.ms rows=1068720
[2025-08-25 07:25:56 UTC] [IViS] [info] [MS] Done: MW-C10_5_MW_chan_-32kms.ms rows=1098366
[2025-08-25 07:25:56 UTC] [IViS] [info] [BLOCK] Done: nchan=1, nbeam=5, nvis_max=1140552 (read 5/5 beams)
Empty sky model¶
sky_model = np.zeros(shape, dtype=np.float32)
Set PB=1 (depends on the response you want: (A*N) or N¶
#pb = np.full(pb.shape,1)
Simulated visibilities (using IViS forward model)¶
#Model visibilities with IVis forward single frequency model
image_processor = Imager3D(I, # visibilities
pb, # array of primary beams
grid, # array of interpolation grids
None, # single dish data in unit of Jy/arcsec^2
None, # beam of single-dish data in radio_beam format
target_header, # header on which to image the data
sky_model, # init array of parameters
0, # maximum number of iterations
0, # hyper-parameter single-dish
False, # impose a positivity constaint
0, # device: 0 is GPU; "cpu" is CPU
0,
beam_workers=1)
model = ClassicIViS3D()
model_vis = image_processor.forward_model(model=model)
[2025-08-25 07:27:26 UTC] [IViS] [info] CUDA unavailable or invalid index; using CPU.
[2025-08-25 07:27:26 UTC] [IViS] [info] CUDA unavailable or invalid index; using CPU.
[2025-08-25 07:27:26 UTC] [IViS] [info] [Initialize Imager3D ]
[2025-08-25 07:27:26 UTC] [IViS] [info] Number of iterations to be performed by the optimizer: 0
[2025-08-25 07:27:26 UTC] [IViS] [warning] lambda_sd = 0 — No short-spacing correction.
[2025-08-25 07:27:26 UTC] [IViS] [info] Optimizer not bounded - Positivity == False
Adding realistic noise from the MeerKAT data beam¶
#Add noise
fact=1 #Scale the noise with this if needed
noise_real = np.random.normal(loc=0.0, scale=I.sigma_I*fact)
noise_imag = np.random.normal(loc=0.0, scale=I.sigma_I*fact)
noise = noise_real + 1j * noise_imag
# Add realistic thermal noise to the model visibilities
I.data_I = model_vis + noise
Imaging using IViS¶
#user parameters
max_its = 20
lambda_sd = 0 #not relevant here
lambda_r = 1 #Control the strength of the Laplacian filtering
cost_device = 0 # 0 for GPU, "cpu" for CPU
optim_device = 0 # 0 for GPU, "cpu" for CPU
positivity = False #Set to False because noise fluctuates around 0
#Initial parameters (zero array)
init_params = np.zeros((1,shape[0],shape[1]), dtype=np.float32)
#create image processor
image_processor = Imager3D(I, # visibilities
pb, # array of primary beams
grid, # array of interpolation grids
sd, # single dish data in unit of Jy/arcsec^2
beam_sd, # beam of single-dish data in radio_beam format
target_header, # header on which to image the data
init_params, # init array of parameters
max_its, # maximum number of iterations
lambda_sd, # hyper-parameter single-dish
positivity, # impose a positivity constaint
cost_device=cost_device,
optim_device=optim_device,
beam_workers=1)
#get image
model = ClassicIViS3D(lambda_r=lambda_r)
noise = image_processor.process(model=model, units="Jy/arcsec^2") #"Jy/arcsec^2" or "K"
[2025-08-25 07:27:28 UTC] [IViS] [info] CUDA unavailable or invalid index; using CPU.
[2025-08-25 07:27:28 UTC] [IViS] [info] CUDA unavailable or invalid index; using CPU.
[2025-08-25 07:27:28 UTC] [IViS] [info] [Initialize Imager3D ]
[2025-08-25 07:27:28 UTC] [IViS] [info] Number of iterations to be performed by the optimizer: 20
[2025-08-25 07:27:28 UTC] [IViS] [warning] lambda_sd = 0 — No short-spacing correction.
[2025-08-25 07:27:28 UTC] [IViS] [info] Optimizer not bounded - Positivity == False
[2025-08-25 07:27:28 UTC] [IViS] [info] Units of output: Jy/arcsec^2.
[2025-08-25 07:27:29 UTC] [IViS] [info] Starting optimisation: PyTorch LBFGS on cpu (unconstrained); cost on cpu
[2025-08-25 07:27:32 UTC] [IViS] [info] [PID 58116] Iter cost: 5.076049e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:27:33 UTC] [IViS] [info] [PID 58116] Iter cost: 5.055916e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:27:35 UTC] [IViS] [info] [PID 58116] Iter cost: 5.030572e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:27:36 UTC] [IViS] [info] [PID 58116] Iter cost: 5.026498e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:27:38 UTC] [IViS] [info] [PID 58116] Iter cost: 5.012884e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:27:40 UTC] [IViS] [info] [PID 58116] Iter cost: 5.008484e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:27:41 UTC] [IViS] [info] [PID 58116] Iter cost: 5.003456e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:27:42 UTC] [IViS] [info] [PID 58116] Iter cost: 4.998165e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:27:44 UTC] [IViS] [info] [PID 58116] Iter cost: 4.995833e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:27:46 UTC] [IViS] [info] [PID 58116] Iter cost: 4.993194e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:27:47 UTC] [IViS] [info] [PID 58116] Iter cost: 4.991972e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:27:49 UTC] [IViS] [info] [PID 58116] Iter cost: 4.990293e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:27:50 UTC] [IViS] [info] [PID 58116] Iter cost: 4.988723e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:27:51 UTC] [IViS] [info] [PID 58116] Iter cost: 4.987264e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:27:53 UTC] [IViS] [info] [PID 58116] Iter cost: 4.986696e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:27:54 UTC] [IViS] [info] [PID 58116] Iter cost: 4.985912e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:27:55 UTC] [IViS] [info] [PID 58116] Iter cost: 4.985510e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:27:57 UTC] [IViS] [info] [PID 58116] Iter cost: 4.984940e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:27:58 UTC] [IViS] [info] [PID 58116] Iter cost: 4.984494e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:27:59 UTC] [IViS] [info] [PID 58116] Iter cost: 4.984227e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:00 UTC] [IViS] [info] [PID 58116] Iter cost: 4.983754e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:00 UTC] [IViS] [info] [Timing] LBFGS (optim_dev=cpu, cost_dev=cpu) took 29.97 s; final loss=5.07605e+06
[2025-08-25 07:28:00 UTC] [IViS] [warning] If you are using ASKAP's convention I = XX + YY (with no 1/2 factor) then multiply the output by 2. Here assuming I = 1/2 (XX + YY).
[2025-08-25 07:28:00 UTC] [IViS] [info] Successful run. Please clap.
#mean pb
filenames = sorted(glob.glob(path_beams+"*.fits"))
n_beams = len(filenames)
pb_all = np.zeros((n_beams,noise.shape[1],noise.shape[2]))
w = dutils.wcs2D(target_header)
for i in tqdm(np.arange(n_beams)):
#open beam cube
hdu_pb = fits.open(filenames[i])
hdr_pb = hdu_pb[0].header
pb2 = hdu_pb[0].data
pb2[pb2 != pb2] = 0.
w_pb = dutils.wcs2D(hdr_pb)
pb2, footprint = reproject_interp((pb2,w_pb.to_header()), w.to_header(), shape)
pb2[pb2 != pb2] = 0.
pb_all += pb2
pb_mean = np.nanmean(pb_all,0)
pb_mean /= np.nanmax(pb_mean)
mask = np.where(pb_mean > 0.05, 1, np.nan)
100%|█████████████████████████████████████████████| 5/5 [00:05<00:00, 1.02s/it]
#PLOT RESULT
fig = plt.figure(figsize=(5, 5))
ax = fig.add_axes([0.1,0.1,0.78,0.8], projection=w)
ax.set_xlabel(r"RA (deg)", fontsize=18.)
ax.set_ylabel(r"DEC (deg)", fontsize=18.)
img = ax.imshow(noise[0]*mask, vmin=-2.5e-5, vmax=2.5e-5, origin="lower")
ax.contour(pb_mean, linestyles="--", levels=[0.2, 0.3], colors=["w","w"])
colorbar_ax = fig.add_axes([0.89, 0.11, 0.02, 0.78])
cbar = fig.colorbar(img, cax=colorbar_ax)
cbar.ax.tick_params(labelsize=14.)
cbar.set_label(r"$I$ (Jy/arcsec$^{2})$", fontsize=18.)
# plt.savefig(pathout + 'ivis_result_cloud_MeerKAT.png', format='png', bbox_inches='tight', pad_inches=0.02, dpi=400)
Compute the P(k)¶
# Get tapper for apodization
tapper = fourier.apodize(0.97, shape)
# Reove the mean
field_zm = noise[0] - np.mean(noise[0])
# Apply apodization function
field_zm_apod = field_zm * tapper
# Compute the 1D P(k)
ks, sps1d_n = fourier.powspec(field_zm_apod, reso=(target_header["CDELT2"]*u.deg).to(u.arcmin).value)
#Plot sps1D
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111)
ax.set_xscale('log')
ax.set_yscale('log')
ax.plot(ks, sps1d_n, "-k")
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel(r"$k$ (arcmin$^{-1}$)", fontsize = 16)
ax.set_ylabel(r"$P(k)$ [(Jy/arcsec$^2$)$^2$]", fontsize = 16)
ax.xaxis.set_major_formatter(FormatStrFormatter('%.3f'))
#plt.savefig('plot/SPS_.png', format='png', bbox_inches='tight', pad_inches=0.02)
Important note¶
It is useful to note that the effective beam response is inherently built into the imaging process. In this case, the sky model contains noise rather than a single point source, but the resulting beam still reflects the non-linear nature of the inversion.
When modeling the power spectrum P(k), the effective beam—shaped by the regularization—acts after the noise, not before. This distinction is important when interpreting the noise properties of the reconstructed image.
This is fundamentally different from a CLEAN approach, where the model is convolved with a chosen restoring beam, and the noise is added back afterward by gridding the residual visibilities.
Compare with P(k) of signal¶
# -------------------
# Read visibilities into VisIData dataclass
# -------------------
reader = CasacoreReader(
prefer_weight_spectrum=False,
keep_autocorr=False,
n_workers=4)
I_s: VisIData = reader.read_blocks_I(
path_ms,
uvmin=0.0,
uvmax=np.inf,
chan_sel=slice(0, 1),
)
# -------------------
# Create Imager3D
# -------------------
image_processor = Imager3D(
vis_data=I_s,
pb=pb,
grid=grid,
sd=sd,
beam_sd=beam_sd,
hdr=target_header,
init_params=init_params,
max_its=max_its,
lambda_sd=lambda_sd,
positivity=positivity,
cost_device=cost_device,
optim_device=optim_device,
beam_workers=1
)
# -------------------
# Choose model
# -------------------
model = ClassicIViS3D(lambda_r=lambda_r, Nw=0)
# get image
result = image_processor.process(model=model, units="Jy/arcsec^2") #"Jy/arcsec^2" or "K"
[2025-08-25 07:28:09 UTC] [IViS] [info] [BLOCK] Loading block from: ./data_tutorials/ivis_data/msl_mw/
[2025-08-25 07:28:09 UTC] [IViS] [info] [BLOCK] Loading 5 beam(s) from: ./data_tutorials/ivis_data/msl_mw/
[2025-08-25 07:28:09 UTC] [IViS] [info] [BLOCK] Parallel read with 4 workers (order-preserving; selective)
[2025-08-25 07:28:09 UTC] [IViS] [info] [MS] Opening: ./data_tutorials/ivis_data/msl_mw/MW-C10_1_MW_chan_-32kms.ms
[2025-08-25 07:28:09 UTC] [IViS] [info] [MS] Opening: ./data_tutorials/ivis_data/msl_mw/MW-C10_4_MW_chan_-32kms.ms
[2025-08-25 07:28:09 UTC] [IViS] [info] [MS] Opening: ./data_tutorials/ivis_data/msl_mw/MW-C10_3_MW_chan_-32kms.ms
[2025-08-25 07:28:09 UTC] [IViS] [info] [MS] Opening: ./data_tutorials/ivis_data/msl_mw/MW-C10_2_MW_chan_-32kms.ms
[2025-08-25 07:28:09 UTC] [IViS] [info] [MS] Done: MW-C10_1_MW_chan_-32kms.ms rows=626313
[2025-08-25 07:28:09 UTC] [IViS] [info] [MS] Opening: ./data_tutorials/ivis_data/msl_mw/MW-C10_5_MW_chan_-32kms.ms
[2025-08-25 07:28:09 UTC] [IViS] [info] [MS] Done: MW-C10_2_MW_chan_-32kms.ms rows=1068720
[2025-08-25 07:28:09 UTC] [IViS] [info] [MS] Done: MW-C10_4_MW_chan_-32kms.ms rows=1140552
[2025-08-25 07:28:09 UTC] [IViS] [info] [MS] Done: MW-C10_3_MW_chan_-32kms.ms rows=1140552
[2025-08-25 07:28:09 UTC] [IViS] [info] [MS] Done: MW-C10_5_MW_chan_-32kms.ms rows=1098366
[2025-08-25 07:28:10 UTC] [IViS] [info] [BLOCK] Done: nchan=1, nbeam=5, nvis_max=1140552 (read 5/5 beams)
[2025-08-25 07:28:10 UTC] [IViS] [info] CUDA unavailable or invalid index; using CPU.
[2025-08-25 07:28:10 UTC] [IViS] [info] CUDA unavailable or invalid index; using CPU.
[2025-08-25 07:28:10 UTC] [IViS] [info] [Initialize Imager3D ]
[2025-08-25 07:28:10 UTC] [IViS] [info] Number of iterations to be performed by the optimizer: 20
[2025-08-25 07:28:10 UTC] [IViS] [warning] lambda_sd = 0 — No short-spacing correction.
[2025-08-25 07:28:10 UTC] [IViS] [info] Optimizer not bounded - Positivity == False
[2025-08-25 07:28:10 UTC] [IViS] [info] Units of output: Jy/arcsec^2.
[2025-08-25 07:28:10 UTC] [IViS] [info] Starting optimisation: PyTorch LBFGS on cpu (unconstrained); cost on cpu
[2025-08-25 07:28:11 UTC] [IViS] [info] [PID 58116] Iter cost: 2.051792e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:12 UTC] [IViS] [info] [PID 58116] Iter cost: 1.943055e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:13 UTC] [IViS] [info] [PID 58116] Iter cost: 1.884155e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:14 UTC] [IViS] [info] [PID 58116] Iter cost: 1.874843e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:16 UTC] [IViS] [info] [PID 58116] Iter cost: 1.868133e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:17 UTC] [IViS] [info] [PID 58116] Iter cost: 1.862304e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:18 UTC] [IViS] [info] [PID 58116] Iter cost: 1.858662e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:19 UTC] [IViS] [info] [PID 58116] Iter cost: 1.855631e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:21 UTC] [IViS] [info] [PID 58116] Iter cost: 1.853403e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:22 UTC] [IViS] [info] [PID 58116] Iter cost: 1.851784e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:23 UTC] [IViS] [info] [PID 58116] Iter cost: 1.850396e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:24 UTC] [IViS] [info] [PID 58116] Iter cost: 1.849463e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:25 UTC] [IViS] [info] [PID 58116] Iter cost: 1.848760e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:27 UTC] [IViS] [info] [PID 58116] Iter cost: 1.848084e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:28 UTC] [IViS] [info] [PID 58116] Iter cost: 1.847788e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:29 UTC] [IViS] [info] [PID 58116] Iter cost: 1.847216e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:30 UTC] [IViS] [info] [PID 58116] Iter cost: 1.847008e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:31 UTC] [IViS] [info] [PID 58116] Iter cost: 1.846706e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:32 UTC] [IViS] [info] [PID 58116] Iter cost: 1.846417e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:33 UTC] [IViS] [info] [PID 58116] Iter cost: 1.846119e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:35 UTC] [IViS] [info] [PID 58116] Iter cost: 1.845987e+06 (optim_dev=cpu, cost_dev=cpu)
[2025-08-25 07:28:35 UTC] [IViS] [info] [Timing] LBFGS (optim_dev=cpu, cost_dev=cpu) took 24.82 s; final loss=2.05179e+06
[2025-08-25 07:28:35 UTC] [IViS] [warning] If you are using ASKAP's convention I = XX + YY (with no 1/2 factor) then multiply the output by 2. Here assuming I = 1/2 (XX + YY).
[2025-08-25 07:28:35 UTC] [IViS] [info] Successful run. Please clap.
#PLOT RESULT
fig = plt.figure(figsize=(5, 5))
ax = fig.add_axes([0.1,0.1,0.78,0.8], projection=w)
ax.set_xlabel(r"RA (deg)", fontsize=18.)
ax.set_ylabel(r"DEC (deg)", fontsize=18.)
img = ax.imshow(result[0]*mask, vmin=-1.e-5, vmax=2.5e-5, origin="lower", cmap="gray_r")
ax.contour(pb_mean, linestyles="--", levels=[0.05, 0.2], colors=["w","w"])
colorbar_ax = fig.add_axes([0.89, 0.11, 0.02, 0.78])
cbar = fig.colorbar(img, cax=colorbar_ax)
cbar.ax.tick_params(labelsize=14.)
cbar.set_label(r"$I$ (Jy/arcsec$^{2})$", fontsize=18.)
# plt.savefig(pathout + 'ivis_result_cloud_MeerKAT.png', format='png', bbox_inches='tight', pad_inches=0.02, dpi=400)
# Get tapper for apodization
tapper = fourier.apodize(0.97, shape)
# Reove the mean
field_zm = result[0] - np.mean(result[0])
# Apply apodization function
field_zm_apod = field_zm * tapper
# Compute the 1D P(k)
ks, sps1d_s = fourier.powspec(field_zm_apod, reso=(target_header["CDELT2"]*u.deg).to(u.arcmin).value)
Rescaling to match amplitude¶
For reasons that I don’t fully understand, the noise realization obtained from synthetizing random noise with SIGMA from the ms doesn’t yield the right scaling. Here this scaling factor is being inferred at ks > kmin_fit (in arcmin).
# Assume ks, sps1d_s, and sps1d_n are already defined
# ---- STEP 1: Define a fitting range (e.g., high-k range) ----
kmin_fit = 0.7 # in arcmin^{-1}, adjust as appropriate
fit_mask = ks > kmin_fit
# ---- STEP 2: Fit scaling factor ----
numerator = np.sum(sps1d_s[fit_mask] * sps1d_n[fit_mask])
denominator = np.sum(sps1d_n[fit_mask] ** 2)
scaling_factor = numerator / denominator
print(f"Best-fit noise scaling factor: {scaling_factor:.4f}")
# ---- STEP 3: Plot ----
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111)
ax.set_xscale('log')
ax.set_yscale('log')
ax.plot(ks, sps1d_s, "-k", label="Signal")
ax.plot(ks, sps1d_n * scaling_factor, "-m", label=f"Noise * {scaling_factor:.2f}")
ax.set_xlabel(r"$k$ (arcmin$^{-1}$)", fontsize=16)
ax.set_ylabel(r"$P(k)$ [(Jy/arcsec$^2$)$^2$]", fontsize=16)
ax.xaxis.set_major_formatter(FormatStrFormatter('%.3f'))
ax.legend()
Best-fit noise scaling factor: 0.3467
<matplotlib.legend.Legend at 0x2a6eb7550>
Here again you can change the regularization strength to see how it both impact the imaged signal and synthetized noise passed through the IViS.