Skip to content

Commit

Permalink
Merge pull request #25 from iniyannatarajan/master
Browse files Browse the repository at this point in the history 
Update noise generation
  • Loading branch information
@iniyannatarajan
iniyannatarajan authored Feb 12, 2024
2 parents 1f7f363 + 37a614a commit 94377db
Show file tree
Hide file tree
Showing 2 changed files with 90 additions and 33 deletions.
18 changes: 4 additions & 14 deletions meqsilhouette/driver/run_meqsilhouette.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -203,7 +203,7 @@ def run_meqsilhouette(config=None):

### TROPOSPHERE COMPONENTS ###
combined_phase_errors = 0 #init for trop combo choice
additive_noises = None

if parameters['trop_enabled']:
info('Tropospheric module is enabled, applying corruptions...')
if parameters['trop_wetonly']:
Expand All @@ -215,16 +215,11 @@ def run_meqsilhouette(config=None):
info('TROPOSPHERE ATTENUATE: attenuating signal using PWV-derived opacity...')
sim_coord.trop_opacity_attenuate()

if parameters['trop_noise']:
info('TROPOSPHERE NOISE: adding sky noise from non-zero PWV...')
additive_noises = sim_coord.trop_add_sky_noise()

if parameters['trop_mean_delay']:
info('TROPOSPHERE DELAY: computing mean delay (time-variability from elevation changes)...')
sim_coord.trop_calc_mean_delays()
combined_phase_errors += sim_coord.phasedelay_alltimes


if parameters['trop_turbulence']:
info('TROPOSPHERE TURBULENCE: computing Kolmogorov turbulence phase errors...')
sim_coord.trop_generate_turbulence_phase_errors()
Expand All @@ -245,9 +240,6 @@ def run_meqsilhouette(config=None):
sim_coord.trop_plots()
info('Generated troposphere plots')

### POPULATE MS WITH SIGMA AND WEIGHT ESTIMATORS ###
sim_coord.add_weights(additive_noises)

### PARALLACTIC ANGLE AND POLARIZATION LEAKAGE ###
if parameters['uvjones_d_on']:
info('Introducing parallactic angle rotation and polarization leakage effects')
Expand All @@ -271,11 +263,9 @@ def run_meqsilhouette(config=None):
info('Generating bandpass plots...')
sim_coord.make_bandpass_plots()

### THERMAL NOISE ###
if parameters['add_thermal_noise']:
info('Adding thermal noise...')
sim_coord.add_receiver_noise()
info('Thermal noise added.')
### Add all noise components and fill in the weight columns
# TODO must do the job of trop_add_sky_noise(), add_weights(additive_noises), add_receiver_noise() here
sim_coord.add_noise(parameters['trop_noise'], parameters['add_thermal_noise'])

### IMAGING, PLOTTING, DATA EXPORT ###
if parameters['make_image']:
Expand Down
105 changes: 86 additions & 19 deletions meqsilhouette/framework/SimCoordinator.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,28 +114,25 @@ def __init__(self, msname, output_column, input_fitsimage, input_fitspol, input_
self.thermal_noise_enabled = thermal_noise_enabled
self.receiver_rms = np.zeros(self.data.shape, dtype='float')

self.compute_weights()

tab.close() # close main MS table


### troposphere information
self.trop_enabled = trop_enabled
if (self.trop_enabled):
self.trop_wetonly = trop_wetonly
self.average_pwv = pwv
self.average_gpress = gpress
self.average_gtemp = gtemp
self.coherence_time = coherence_time
self.fixdelay_max_picosec = fixdelay_max_picosec
self.elevation_tropshape = np.expand_dims(np.swapaxes(self.elevation, 0, 1), 1) # reshaped for troposphere operations
self.opacity, self.sky_temp = self.trop_return_opacity_sky_temp()
self.transmission = np.exp(-1*self.opacity)

# Set some optional arrays to None. These will be filled later depending upon the user request.
self.transmission_matrix = None
self.turb_phase_errors = None
self.delay_alltimes = None
self.sky_noise = None
self.trop_wetonly = trop_wetonly
self.average_pwv = pwv
self.average_gpress = gpress
self.average_gtemp = gtemp
self.coherence_time = coherence_time
self.fixdelay_max_picosec = fixdelay_max_picosec
self.elevation_tropshape = np.expand_dims(np.swapaxes(self.elevation, 0, 1), 1) # reshaped for troposphere operations
self.opacity, self.sky_temp = self.trop_return_opacity_sky_temp()
self.transmission = np.exp(-1*self.opacity)

# Set some optional arrays to None. These will be filled later depending upon the user request.
self.transmission_matrix = None
self.turb_phase_errors = None
self.delay_alltimes = None
self.sky_noise = None

### bandpass information
self.bandpass_table = bandpass_table
Expand Down Expand Up @@ -1091,6 +1088,76 @@ def add_gjones_manual(self):

self.save_data()

##################################
# Add noise components
def add_noise(self, tropnoise, thermalnoise):
""" Add sky and receiver noise components to the visibilities and populate weight columns"""

# compute receiver rms noise
for a0 in range(self.Nant):
for a1 in range(self.Nant):
if a1 > a0:
self.receiver_rms[self.baseline_dict[(a0,a1)]] = (1/self.corr_eff) * np.sqrt(self.SEFD[a0] * \
self.SEFD[a1] / float(2 * self.tint * self.chan_width))

# compute sky sigma estimator (i.e. sky rms noise) and realise sky noise
if tropnoise:
sefd_matrix = 2 * Boltzmann / self.dish_area * (1e26*self.sky_temp * (1. - np.exp(-1.0 * self.opacity / np.sin(self.elevation_tropshape))))
self.sky_noise = np.zeros(self.data.shape, dtype='complex')
sky_sigma_estimator = np.zeros(self.data.shape)

for a0 in range(self.Nant):
for a1 in range(self.Nant):
if a1 > a0:
rms = (1/self.corr_eff) * np.sqrt(sefd_matrix[:, :, a0] * sefd_matrix[:, :, a1] / (float(2 * self.tint * self.chan_width)))
self.temp_rms = rms
rms = np.expand_dims(rms, 2)
rms = rms * np.ones((1, 1, 4))
self.sky_noise[self.baseline_dict[(a0, a1)]] = self.rng_atm.normal(0.0, rms) + 1j * self.rng_atm.normal(0.0, rms)
sky_sigma_estimator[self.baseline_dict[(a0, a1)]] = rms
np.save(II('$OUTDIR')+'/atm_output/sky_noise_timestamp_%d'%(self.timestamp), self.sky_noise)

# add sky noise rms to receiver rms if tropnoise is set
try:
for tind in range(self.nchunks):
self.receiver_rms[tind*self.chunksize:(tind+1)*self.chunksize] = np.sqrt(np.power(self.receiver_rms[tind*self.chunksize:(tind+1)*self.chunksize], 2) + \
np.power(sky_sigma_estimator[tind*self.chunksize:(tind+1)*self.chunksize], 2))
except MemoryError:
abort("Arrays too large to be held in memory. Aborting execution.")

# use the receiver rms to realise thermal noise
if thermalnoise:
info('Realise thermal noise from receiver rms...')
self.thermal_noise = np.zeros(self.data.shape, dtype='complex')
size = (self.time_unique.shape[0], self.chan_freq.shape[0], 4)
for a0 in range(self.Nant):
for a1 in range(self.Nant):
if a1 > a0:
rms = self.receiver_rms[self.baseline_dict[(a0,a1)]]
self.thermal_noise[self.baseline_dict[(a0, a1)]] = self.rng_predict.normal(0.0, rms, size=size) + 1j * self.rng_predict.normal(0.0, rms, size=size)

np.save(II('$OUTDIR')+'/receiver_noise_timestamp_%d'%(self.timestamp), self.thermal_noise)

# add both the noise components to the data
try:
info('Applying thermal noise to data...')
for tind in range(self.nchunks):
self.data[tind*self.chunksize:(tind+1)*self.chunksize] += self.thermal_noise[tind*self.chunksize:(tind+1)*self.chunksize] + \
self.sky_noise[tind*self.chunksize:(tind+1)*self.chunksize]
self.save_data()
except MemoryError:
abort("Arrays too large to be held in memory. Aborting execution.")

# populate MS weight and sigma columns using the receiver rms (with or without sky noise)
tab = pt.table(self.msname, readonly=False,ack=False)
tab.putcol("SIGMA", self.receiver_rms[:,0,:])
if 'SIGMA_SPECTRUM' in tab.colnames():
tab.putcol("SIGMA_SPECTRUM", self.receiver_rms)
tab.putcol("WEIGHT", 1/self.receiver_rms[:,0,:]**2)
if 'WEIGHT_SPECTRUM' in tab.colnames():
tab.putcol("WEIGHT_SPECTRUM", 1/self.receiver_rms**2)
tab.close()

#############################
##### General MS plots #####
#############################
Expand Down

0 comments on commit 94377db

Please sign in to comment.