Skip to content

Commit a7c434b

Browse files
timcornwelltimcornwell
committed
Best arl-demo so far
1 parent a7356f2 commit a7c434b

22 files changed

Lines changed: 495 additions & 198 deletions

‎processing_components/calibration/calibration.py‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -46,8 +46,7 @@ def solve_gaintable(vis: BlockVisibility, modelvis: BlockVisibility = None, gt=N
4646
assert isinstance(vis, BlockVisibility), vis
4747
if modelvis is not None:
4848
assert isinstance(modelvis, BlockVisibility), modelvis
49-
50-
assert numpy.max(numpy.abs(modelvis.vis)) > 0.0, "Model visibility is zero"
49+
assert numpy.max(numpy.abs(modelvis.vis)) > 0.0, "Model visibility is zero"
5150

5251
if phase_only:
5352
log.debug('solve_gaintable: Solving for phase only')

‎processing_components/calibration/calibration_control.py‎

Lines changed: 80 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -46,11 +46,45 @@ def create_calibration_controls(**kwargs):
4646
return controls
4747

4848

49-
def calibrate_function(vis, model_vis, calibration_context='T', controls=None, iteration=0, tol=1e-6, **kwargs):
49+
def apply_calibration_function(vis, gaintables, calibration_context='T', controls=None, iteration=0, tol=1e-6,
50+
**kwargs):
5051
""" Calibrate using algorithm specified by calibration_context
5152
5253
The context string can denote a sequence of calibrations e.g. TGB with different timescales.
5354
55+
:param vis:
56+
:param model_vis:
57+
:param calibration_context: calibration contexts in order of correction e.g. 'TGB'
58+
:param control: controls dictionary, modified as necessary
59+
:param iteration: Iteration number to be compared to the 'first_selfcal' field.
60+
:param kwargs:
61+
:return: Calibrated data_models, dict(gaintables)
62+
"""
63+
64+
if controls is None:
65+
controls = create_calibration_controls(**kwargs)
66+
67+
isVis = isinstance(vis, Visibility)
68+
if isVis:
69+
avis = convert_visibility_to_blockvisibility(vis)
70+
else:
71+
avis = vis
72+
73+
for c in calibration_context:
74+
if iteration >= controls[c]['first_selfcal']:
75+
avis = apply_gaintable(avis, gaintables[c], inverse=True, timeslice=controls[c]['timeslice'])
76+
77+
if isVis:
78+
return convert_blockvisibility_to_visibility(avis)
79+
else:
80+
return avis
81+
82+
83+
def calibrate_function(vis, model_vis, calibration_context='T', controls=None, iteration=0, tol=1e-6, **kwargs):
84+
""" Calibrate using algorithm specified by calibration_context
85+
86+
The context string can denote a sequence of calibrations e.g. TGB with different timescales.
87+
5488
:param vis:
5589
:param model_vis:
5690
:param calibration_context: calibration contexts in order of correction e.g. 'TGB'
@@ -94,3 +128,48 @@ def calibrate_function(vis, model_vis, calibration_context='T', controls=None, i
94128
return convert_blockvisibility_to_visibility(avis), gaintables
95129
else:
96130
return avis, gaintables
131+
132+
133+
def solve_calibrate_function(vis, model_vis, calibration_context='T', controls=None, iteration=0, tol=1e-6, **kwargs):
134+
""" Calibrate using algorithm specified by calibration_context
135+
136+
The context string can denote a sequence of calibrations e.g. TGB with different timescales.
137+
138+
:param vis:
139+
:param model_vis:
140+
:param calibration_context: calibration contexts in order of correction e.g. 'TGB'
141+
:param control: controls dictionary, modified as necessary
142+
:param iteration: Iteration number to be compared to the 'first_selfcal' field.
143+
:param kwargs:
144+
:return: Calibrated data_models, dict(gaintables)
145+
"""
146+
gaintables = {}
147+
148+
if controls is None:
149+
controls = create_calibration_controls(**kwargs)
150+
151+
isVis = isinstance(vis, Visibility)
152+
if isVis:
153+
avis = convert_visibility_to_blockvisibility(vis)
154+
else:
155+
avis = vis
156+
157+
isMVis = isinstance(model_vis, Visibility)
158+
if isMVis:
159+
amvis = convert_visibility_to_blockvisibility(model_vis)
160+
else:
161+
amvis = model_vis
162+
163+
for c in calibration_context:
164+
if iteration >= controls[c]['first_selfcal']:
165+
gaintables[c] = \
166+
create_gaintable_from_blockvisibility(avis, timeslice=controls[c]['timeslice'])
167+
gaintables[c] = solve_gaintable(avis, amvis,
168+
timeslice=controls[c]['timeslice'],
169+
phase_only=controls[c]['phase_only'],
170+
crosspol=controls[c]['shape'] == 'matrix',
171+
tol=tol)
172+
log.debug(qa_gaintable(gaintables[c], context='Jones matrix %s, iteration %d' % (c, iteration)))
173+
else:
174+
log.debug('calibrate_function: Jones matrix %s not solved, iteration %d' % (c, iteration))
175+
return gaintables

‎processing_components/calibration/operations.py‎

Lines changed: 13 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -130,18 +130,21 @@ def apply_gaintable(vis: BlockVisibility, gt: GainTable, inverse=False, vis_slic
130130
original = vis.vis[rows]
131131
applied = copy.deepcopy(original)
132132
for time in range(ntimes):
133-
for a1 in range(vis.nants - 1):
134-
for a2 in range(a1 + 1, vis.nants):
135-
if is_scalar:
136-
smueller = gain[time, a1, :, 0] * numpy.conjugate(gain[time, a2, :, 0])
137-
if inverse:
138-
if numpy.abs(smueller).all() > 0.0:
139-
applied[time, a2, a1, :, 0][..., numpy.newaxis] = \
140-
original[time, a2, a1, :, 0][..., numpy.newaxis] / smueller
141-
else:
133+
if is_scalar:
134+
if inverse:
135+
lgain = numpy.ones_like(gain)
136+
lgain[numpy.abs(gain) > 0.0] = 1.0 / gain[numpy.abs(gain) > 0.0]
137+
else:
138+
lgain = gain
139+
for a1 in range(vis.nants - 1):
140+
for a2 in range(a1 + 1, vis.nants):
141+
if is_scalar:
142+
smueller = lgain[time, a1, :, 0] * numpy.conjugate(lgain[time, a2, :, 0])
142143
applied[time, a2, a1, :, 0][..., numpy.newaxis] = \
143144
original[time, a2, a1, :, 0][..., numpy.newaxis] * smueller
144-
else:
145+
else:
146+
for a1 in range(vis.nants - 1):
147+
for a2 in range(a1 + 1, vis.nants):
145148
for chan in range(nchan):
146149
mueller = numpy.kron(gain[time, a1, chan, :, :],
147150
numpy.conjugate(gain[time, a2, chan, :, :]))

‎processing_components/image/deconvolution.py‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -59,7 +59,7 @@ def deconvolve_cube(dirty: Image, psf: Image, prefix='', **kwargs) -> (Image, Im
5959
comp, residual = deconvolve_cube(dirty, psf, niter=1000, gain=0.7, algorithm='msclean',
6060
scales=[0, 3, 10, 30], threshold=0.01)
6161
62-
For the MFS clean, the psf must have number of channels >= 2 * nmoments
62+
For the MFS clean, the psf must have number of channels >= 2 * nmoment
6363
6464
:param dirty: Image dirty image
6565
:param psf: Image Point Spread Function
@@ -69,7 +69,7 @@ def deconvolve_cube(dirty: Image, psf: Image, prefix='', **kwargs) -> (Image, Im
6969
:param threshold: Clean threshold (0.0)
7070
:param fractional_threshold: Fractional threshold (0.01)
7171
:param scales: Scales (in pixels) for multiscale ([0, 3, 10, 30])
72-
:param nmoments: Number of frequency moments (default 3)
72+
:param nmoment: Number of frequency moments (default 3)
7373
:param findpeak: Method of finding peak in mfsclean: 'Algorithm1'|'ASKAPSoft'|'CASA'|'ARL', Default is ARL.
7474
:return: componentimage, residual
7575
@@ -137,12 +137,12 @@ def deconvolve_cube(dirty: Image, psf: Image, prefix='', **kwargs) -> (Image, Im
137137

138138
log.info("deconvolve_cube %s: Multi-scale multi-frequency clean of each polarisation separately"
139139
% prefix)
140-
nmoments = get_parameter(kwargs, "nmoments", 3)
141-
assert nmoments > 0, "Number of frequency moments must be greater than zero"
140+
nmoment = get_parameter(kwargs, "nmoment", 3)
141+
assert nmoment > 0, "Number of frequency moments must be greater than zero"
142142
nchan = dirty.shape[0]
143-
assert nchan > 2 * nmoments, "Require nchan %d > 2 * nmoments %d" % (nchan, 2 * nmoments)
144-
dirty_taylor = calculate_image_frequency_moments(dirty, nmoments=nmoments)
145-
psf_taylor = calculate_image_frequency_moments(psf, nmoments=2 * nmoments)
143+
assert nchan > 2 * nmoment, "Require nchan %d > 2 * nmoment %d" % (nchan, 2 * nmoment)
144+
dirty_taylor = calculate_image_frequency_moments(dirty, nmoment=nmoment)
145+
psf_taylor = calculate_image_frequency_moments(psf, nmoment=2 * nmoment)
146146
psf_peak = numpy.max(psf_taylor.data)
147147
dirty_taylor.data /= psf_peak
148148
psf_taylor.data /= psf_peak

‎processing_components/image/operations.py‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -197,7 +197,7 @@ def smooth_image(model: Image, width=1.0):
197197
return cmodel
198198

199199

200-
def calculate_image_frequency_moments(im: Image, reference_frequency=None, nmoments=3) -> Image:
200+
def calculate_image_frequency_moments(im: Image, reference_frequency=None, nmoment=3) -> Image:
201201
"""Calculate frequency weighted moments
202202
203203
Weights are ((freq-reference_frequency)/reference_frequency)**moment
@@ -206,28 +206,28 @@ def calculate_image_frequency_moments(im: Image, reference_frequency=None, nmome
206206
207207
For example, to find the moments and then reconstruct from just the moments::
208208
209-
moment_cube = calculate_image_frequency_moments(model_multichannel, nmoments=5)
209+
moment_cube = calculate_image_frequency_moments(model_multichannel, nmoment=5)
210210
reconstructed_cube = calculate_image_from_frequency_moments(model_multichannel, moment_cube)
211211
212212
:param im: Image cube
213213
:param reference_frequency: Reference frequency (default None uses average)
214-
:param nmoments: Number of moments to calculate
214+
:param nmoment: Number of moments to calculate
215215
:return: Moments image
216216
"""
217217
assert isinstance(im, Image)
218218
nchan, npol, ny, nx = im.shape
219219
channels = numpy.arange(nchan)
220220
freq = im.wcs.sub(['spectral']).wcs_pix2world(channels, 0)[0]
221221

222-
assert nmoments <= nchan, "Number of moments %d cannot exceed the number of channels %d" % (nmoments, nchan)
222+
assert nmoment <= nchan, "Number of moments %d cannot exceed the number of channels %d" % (nmoment, nchan)
223223

224224
if reference_frequency is None:
225225
reference_frequency = numpy.average(freq)
226226
log.debug("calculate_image_frequency_moments: Reference frequency = %.3f (MHz)" % (reference_frequency/1e6))
227227

228-
moment_data = numpy.zeros([nmoments, npol, ny, nx])
228+
moment_data = numpy.zeros([nmoment, npol, ny, nx])
229229

230-
for moment in range(nmoments):
230+
for moment in range(nmoment):
231231
for chan in range(nchan):
232232
weight = numpy.power((freq[chan] - reference_frequency) / reference_frequency, moment)
233233
moment_data[moment, ...] += im.data[chan, ...] * weight
@@ -251,7 +251,7 @@ def calculate_image_from_frequency_moments(im: Image, moment_image: Image, refer
251251
252252
For example, to find the moments and then reconstruct from just the moments::
253253
254-
moment_cube = calculate_image_frequency_moments(model_multichannel, nmoments=5)
254+
moment_cube = calculate_image_frequency_moments(model_multichannel, nmoment=5)
255255
reconstructed_cube = calculate_image_from_frequency_moments(model_multichannel, moment_cube)
256256
257257
@@ -262,7 +262,7 @@ def calculate_image_from_frequency_moments(im: Image, moment_image: Image, refer
262262
"""
263263
assert isinstance(im, Image)
264264
nchan, npol, ny, nx = im.shape
265-
nmoments, mnpol, mny, mnx = moment_image.shape
265+
nmoment, mnpol, mny, mnx = moment_image.shape
266266

267267
assert npol == mnpol
268268
assert ny == mny
@@ -281,7 +281,7 @@ def calculate_image_from_frequency_moments(im: Image, moment_image: Image, refer
281281

282282
newim.data[...] = 0.0
283283

284-
for moment in range(nmoments):
284+
for moment in range(nmoment):
285285
for chan in range(nchan):
286286
weight = numpy.power((freq[chan] - reference_frequency) / reference_frequency, moment)
287287
newim.data[chan, ...] += moment_image.data[moment, ...] * weight

‎processing_components/simulation/testing_support.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -836,13 +836,13 @@ def insert_unittest_errors(vt, seed=180555, calibration_context="TGB", amp_error
836836
phase_errors = {'T': 1.0, 'G': 0.1, 'B': 0.01}
837837

838838
for c in calibration_context:
839-
gaintable = \
840-
create_gaintable_from_blockvisibility(vt, timeslice=controls[c]['timeslice'])
839+
gaintable = create_gaintable_from_blockvisibility(vt, timeslice=controls[c]['timeslice'])
841840
gaintable = simulate_gaintable(gaintable,
842841
timeslice=controls[c]['timeslice'],
843842
phase_only=controls[c]['phase_only'],
844843
crosspol=controls[c]['shape'] == 'matrix',
845844
phase_error=phase_errors[c], amplitude_error=amp_errors[c])
845+
846846
vt = apply_gaintable(vt, gaintable, inverse=True, timeslice=controls[c]['timeslice'])
847847

848848
return vt

‎processing_components/visibility/gather_scatter.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -141,6 +141,7 @@ def visibility_gather_channel(vis_list: List[BlockVisibility], vis: BlockVisibil
141141
for v in vis_list:
142142
assert len(v.frequency) == 1
143143
assert len(v.channel_bandwidth) == 1
144+
144145
vis = BlockVisibility(data=None,
145146
frequency=numpy.array([v.frequency[0] for v in vis_list]),
146147
channel_bandwidth=numpy.array([v.channel_bandwidth[0] for v in vis_list]),

0 commit comments

Comments
 (0)