Skip to content

Commit bded1b6

Browse files
ChrisHadChrisHad
authored
Merge pull request #18 from SKA-ScienceDataProcessor/master
Merge main repo
2 parents e4d31cf + a289c3a commit bded1b6

19 files changed

Lines changed: 331 additions & 267 deletions

File tree

‎Makefile‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -86,7 +86,7 @@ code-lint:
8686
code-analysis: code-flake code-lint
8787

8888
examples: inplace
89-
$(MAKE) -C libs/notebooks
89+
$(MAKE) -C processing_library/notebooks
9090
$(MAKE) -C processing_components/notebooks
9191
$(MAKE) -C workflows/notebooks
9292

‎processing_components/calibration/calibration_control.py‎

Lines changed: 97 additions & 60 deletions
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@
1212

1313
import numpy
1414

15-
from data_models.memory_data_models import Visibility
15+
from data_models.memory_data_models import Visibility, BlockVisibility
1616

1717
from ..calibration.operations import create_gaintable_from_blockvisibility, apply_gaintable, qa_gaintable
1818
from ..calibration.calibration import solve_gaintable
@@ -65,21 +65,33 @@ def apply_calibration_function(vis, gaintables, calibration_context='T', control
6565

6666
if controls is None:
6767
controls = create_calibration_controls(**kwargs)
68-
69-
isVis = isinstance(vis, Visibility)
70-
if isVis:
71-
avis = convert_visibility_to_blockvisibility(vis)
72-
else:
73-
avis = vis
74-
68+
69+
# Check to see if changes are required
70+
changes = False
7571
for c in calibration_context:
7672
if iteration >= controls[c]['first_selfcal']:
77-
avis = apply_gaintable(avis, gaintables[c], inverse=True, timeslice=controls[c]['timeslice'])
73+
changes = True
74+
75+
if changes:
76+
77+
isVis = isinstance(vis, Visibility)
78+
if isVis:
79+
avis = convert_visibility_to_blockvisibility(vis)
80+
else:
81+
avis = vis
82+
83+
assert isinstance(avis, BlockVisibility), avis
84+
85+
for c in calibration_context:
86+
if iteration >= controls[c]['first_selfcal']:
87+
avis = apply_gaintable(avis, gaintables[c], inverse=True, timeslice=controls[c]['timeslice'])
7888

79-
if isVis:
80-
return convert_blockvisibility_to_visibility(avis)
89+
if isVis:
90+
return convert_blockvisibility_to_visibility(avis)
91+
else:
92+
return avis
8193
else:
82-
return avis
94+
return vis
8395

8496

8597
def calibrate_function(vis, model_vis, calibration_context='T', controls=None, iteration=0, tol=1e-8, **kwargs):
@@ -100,36 +112,50 @@ def calibrate_function(vis, model_vis, calibration_context='T', controls=None, i
100112
if controls is None:
101113
controls = create_calibration_controls(**kwargs)
102114

103-
isVis = isinstance(vis, Visibility)
104-
if isVis:
105-
avis = convert_visibility_to_blockvisibility(vis)
106-
else:
107-
avis = vis
108-
109-
isMVis = isinstance(model_vis, Visibility)
110-
if isMVis:
111-
amvis = convert_visibility_to_blockvisibility(model_vis)
112-
else:
113-
amvis = model_vis
114-
115+
# Check to see if changes are required
116+
changes = False
115117
for c in calibration_context:
116118
if iteration >= controls[c]['first_selfcal']:
117-
gaintables[c] = \
118-
create_gaintable_from_blockvisibility(avis, timeslice=controls[c]['timeslice'])
119-
gaintables[c] = solve_gaintable(avis, amvis,
120-
timeslice=controls[c]['timeslice'],
121-
phase_only=controls[c]['phase_only'],
122-
crosspol=controls[c]['shape'] == 'matrix',
123-
tol=tol)
124-
log.debug('calibrate_function: Jones matrix %s, iteration %d' % (c, iteration))
125-
log.debug(qa_gaintable(gaintables[c], context='Jones matrix %s, iteration %d' % (c, iteration)))
126-
avis = apply_gaintable(avis, gaintables[c], inverse=True, timeslice=controls[c]['timeslice'])
119+
changes = True
120+
121+
if changes:
122+
123+
isVis = isinstance(vis, Visibility)
124+
if isVis:
125+
avis = convert_visibility_to_blockvisibility(vis)
127126
else:
128-
log.debug('calibrate_function: Jones matrix %s not solved, iteration %d' % (c, iteration))
129-
if isVis:
130-
return convert_blockvisibility_to_visibility(avis), gaintables
127+
avis = vis
128+
129+
isMVis = isinstance(model_vis, Visibility)
130+
if isMVis:
131+
amvis = convert_visibility_to_blockvisibility(model_vis)
132+
else:
133+
amvis = model_vis
134+
135+
assert isinstance(avis, BlockVisibility), avis
136+
137+
for c in calibration_context:
138+
if iteration >= controls[c]['first_selfcal']:
139+
changed = True
140+
gaintables[c] = \
141+
create_gaintable_from_blockvisibility(avis, timeslice=controls[c]['timeslice'])
142+
gaintables[c] = solve_gaintable(avis, amvis,
143+
timeslice=controls[c]['timeslice'],
144+
phase_only=controls[c]['phase_only'],
145+
crosspol=controls[c]['shape'] == 'matrix',
146+
tol=tol)
147+
log.debug('calibrate_function: Jones matrix %s, iteration %d' % (c, iteration))
148+
log.debug(qa_gaintable(gaintables[c], context='Jones matrix %s, iteration %d' % (c, iteration)))
149+
avis = apply_gaintable(avis, gaintables[c], inverse=True, timeslice=controls[c]['timeslice'])
150+
else:
151+
log.debug('calibrate_function: Jones matrix %s not solved, iteration %d' % (c, iteration))
152+
153+
if isVis:
154+
return convert_blockvisibility_to_visibility(avis), gaintables
155+
else:
156+
return avis, gaintables
131157
else:
132-
return avis, gaintables
158+
return vis, gaintables
133159

134160

135161
def solve_calibrate_function(vis, model_vis, calibration_context='T', controls=None, iteration=0, tol=1e-6, **kwargs):
@@ -149,29 +175,40 @@ def solve_calibrate_function(vis, model_vis, calibration_context='T', controls=N
149175

150176
if controls is None:
151177
controls = create_calibration_controls(**kwargs)
152-
153-
isVis = isinstance(vis, Visibility)
154-
if isVis:
155-
avis = convert_visibility_to_blockvisibility(vis)
156-
else:
157-
avis = vis
158-
159-
isMVis = isinstance(model_vis, Visibility)
160-
if isMVis:
161-
amvis = convert_visibility_to_blockvisibility(model_vis)
162-
else:
163-
amvis = model_vis
164-
178+
179+
# Check to see if changes are required
180+
changes = False
165181
for c in calibration_context:
166182
if iteration >= controls[c]['first_selfcal']:
167-
gaintables[c] = \
168-
create_gaintable_from_blockvisibility(avis, timeslice=controls[c]['timeslice'])
169-
gaintables[c] = solve_gaintable(avis, amvis,
170-
timeslice=controls[c]['timeslice'],
171-
phase_only=controls[c]['phase_only'],
172-
crosspol=controls[c]['shape'] == 'matrix',
173-
tol=tol)
174-
log.debug(qa_gaintable(gaintables[c], context='Jones matrix %s, iteration %d' % (c, iteration)))
183+
changes = True
184+
185+
if changes:
186+
187+
isVis = isinstance(vis, Visibility)
188+
if isVis:
189+
avis = convert_visibility_to_blockvisibility(vis)
175190
else:
176-
log.debug('calibrate_function: Jones matrix %s not solved, iteration %d' % (c, iteration))
191+
avis = vis
192+
193+
isMVis = isinstance(model_vis, Visibility)
194+
if isMVis:
195+
amvis = convert_visibility_to_blockvisibility(model_vis)
196+
else:
197+
amvis = model_vis
198+
199+
assert isinstance(avis, BlockVisibility), avis
200+
201+
for c in calibration_context:
202+
if iteration >= controls[c]['first_selfcal']:
203+
gaintables[c] = \
204+
create_gaintable_from_blockvisibility(avis, timeslice=controls[c]['timeslice'])
205+
gaintables[c] = solve_gaintable(avis, amvis,
206+
timeslice=controls[c]['timeslice'],
207+
phase_only=controls[c]['phase_only'],
208+
crosspol=controls[c]['shape'] == 'matrix',
209+
tol=tol)
210+
log.debug(qa_gaintable(gaintables[c], context='Jones matrix %s, iteration %d' % (c, iteration)))
211+
else:
212+
log.debug('calibrate_function: Jones matrix %s not solved, iteration %d' % (c, iteration))
213+
177214
return gaintables

‎processing_components/simulation/testing_support.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -826,7 +826,6 @@ def insert_unittest_errors(vt, seed=180555, calibration_context="TG", amp_errors
826826
:param amp_errors: e.g. {'T': 0.0, 'G': 0.01, 'B': 0.01}
827827
:return:
828828
"""
829-
numpy.random.seed(seed)
830829
controls = create_calibration_controls()
831830

832831
if amp_errors is None:
@@ -841,7 +840,8 @@ def insert_unittest_errors(vt, seed=180555, calibration_context="TG", amp_errors
841840
timeslice=controls[c]['timeslice'],
842841
phase_only=controls[c]['phase_only'],
843842
crosspol=controls[c]['shape'] == 'matrix',
844-
phase_error=phase_errors[c], amplitude_error=amp_errors[c])
843+
phase_error=phase_errors[c], amplitude_error=amp_errors[c],
844+
seed=seed)
845845

846846
vt = apply_gaintable(vt, gaintable, inverse=True, timeslice=controls[c]['timeslice'])
847847

‎processing_library/calibration/solvers.py‎

Lines changed: 18 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -56,7 +56,8 @@ def solve_from_X(gt: GainTable, x: numpy.ndarray, xwt: numpy.ndarray, chunk, cro
5656
return gt
5757

5858

59-
def solve_antenna_gains_itsubs_scalar(gain, gwt, x, xwt, niter=30, tol=1e-8, phase_only=True, refant=0):
59+
def solve_antenna_gains_itsubs_scalar(gain, gwt, x, xwt, niter=30, tol=1e-8, phase_only=True, refant=0,
60+
damping=0.5):
6061
"""Solve for the antenna gains
6162
6263
x(antenna2, antenna1) = gain(antenna1) conj(gain(antenna2))
@@ -93,11 +94,11 @@ def solve_antenna_gains_itsubs_scalar(gain, gwt, x, xwt, niter=30, tol=1e-8, pha
9394
gain[mask] = gain[mask] / numpy.abs(gain[mask])
9495
angles = numpy.angle(gain)
9596
gain *= numpy.exp(-1j * angles)[refant, ...]
96-
gain = 0.5 * (gain + gainLast)
97+
gain = (1.0 - damping) * gain + damping * gainLast
9798
change = numpy.max(numpy.abs(gain - gainLast))
9899
if change < tol:
99100
return gain, gwt, solution_residual_scalar(gain, x, xwt)
100-
101+
101102
return gain, gwt, solution_residual_scalar(gain, x, xwt)
102103

103104

@@ -113,7 +114,7 @@ def gain_substitution_scalar(gain, x, xwt):
113114
top = numpy.sum(x[:, ant1, :, 0, 0] * gain[:, :, 0, 0] * xwt[:, ant1, :, 0, 0], axis=0)
114115
bot = numpy.sum((gain[:, :, 0, 0] * numpy.conjugate(gain[:, :, 0, 0]) *
115116
xwt[:, ant1, :, 0, 0]).real, axis=0)
116-
117+
117118
if bot.all() > 0.0:
118119
newgain[ant1, :, 0, 0] = top / bot
119120
gwt[ant1, :, 0, 0] = bot
@@ -304,22 +305,24 @@ def solution_residual_scalar(gain, x, xwt):
304305
:return: residual[...]
305306
"""
306307

307-
nants, nchan, nrec, _ = gain.shape
308-
x = x.reshape(nants, nants, nchan, nrec, nrec)
308+
nant, nchan, nrec, _ = gain.shape
309+
x = x.reshape(nant, nant, nchan, nrec, nrec)
309310

310-
xwt = xwt.reshape(nants, nants, nchan, nrec, nrec)
311+
xwt = xwt.reshape(nant, nant, nchan, nrec, nrec)
311312

312313
residual = numpy.zeros([nchan, nrec, nrec])
313314
sumwt = numpy.zeros([nchan, nrec, nrec])
314315

315-
for ant1 in range(nants):
316-
for ant2 in range(nants):
317-
for chan in range(nchan):
318-
error = x[ant2, ant1, chan, 0, 0] - \
319-
gain[ant1, chan, 0, 0] * numpy.conjugate(gain[ant2, chan, 0, 0])
320-
residual += (error * xwt[ant2, ant1, chan, 0, 0] * numpy.conjugate(error)).real
321-
sumwt += xwt[ant2, ant1, chan, 0, 0]
322-
316+
for chan in range(nchan):
317+
lgain = gain[:, chan, 0, 0]
318+
clgain = numpy.conjugate(lgain)
319+
smueller = numpy.ma.outer(clgain, lgain).reshape([nant, nant])
320+
error = x[:, :, chan, 0, 0] - smueller
321+
for i in range(nant):
322+
error[i,i]=0.0
323+
residual += numpy.sum(error * xwt[:, :, chan, 0, 0] * numpy.conjugate(error)).real
324+
sumwt += numpy.sum(xwt[:, :, chan, 0, 0])
325+
323326
residual[sumwt > 0.0] = numpy.sqrt(residual[sumwt > 0.0] / sumwt[sumwt > 0.0])
324327
residual[sumwt <= 0.0] = 0.0
325328

‎processing_library/util/coordinate_support.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -185,7 +185,7 @@ def lmn_to_skycoord(lmn, phasecentre: SkyCoord):
185185
# Convert l,m,n to SkyCoord convention, also enforce celestial sphere
186186
n = numpy.sqrt(1 - lmn[0] ** 2 - lmn[1] ** 2) - 1.0
187187
dc = n + 1, lmn[0], lmn[1]
188-
target = SkyCoord(x=dc[0], y=dc[1], z=dc[2], representation='cartesian', frame=phasecentre.skyoffset_frame())
188+
target = SkyCoord(x=dc[0], y=dc[1], z=dc[2], representation_type='cartesian', frame=phasecentre.skyoffset_frame())
189189
return target.transform_to(phasecentre.frame)
190190

191191

‎setup.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -42,7 +42,7 @@
4242
author_email='realtimcornwell@gmail.com',
4343
url='https://github.com/SKA-ScienceDataProcessor/algorithm-reference-library',
4444
license='Apache License Version 2.0',
45-
packages=['data_models', 'processing_library', 'processing_components', 'workflows', 'examples'],
45+
packages=['data_models', 'processing_library', 'processing_components', 'workflows', 'wrappers', 'examples'],
4646
test_suite="tests",
4747
tests_require=['pytest'],
4848
ext_modules = [libarlffi]

0 commit comments

Comments
 (0)