Skip to content

Commit 419d650

Browse files
achaelachael
committed
added apply_singepol to selfcal
1 parent 7fb8967 commit 419d650

1 file changed

Lines changed: 78 additions & 23 deletions

File tree

‎ehtim/calibrating/self_cal.py‎

Lines changed: 78 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -50,11 +50,12 @@
5050
###################################################################################################
5151

5252

53-
def self_cal(obs, im, sites=[], method="both", pol='I', minimizer_method='BFGS',
53+
def self_cal(obs, im, sites=[], pol='I', apply_singlepol=False, method="both",
54+
minimizer_method='BFGS',
5455
pad_amp=0., gain_tol=.2, solution_interval=0.0, scan_solutions=False,
5556
ttype='direct', fft_pad_factor=2, caltable=False,
5657
debias=True, apply_dterms=False,
57-
copy_closure_tables=True,
58+
copy_closure_tables=False,
5859
processes=-1, show_solution=False, msgtype='bar',
5960
use_grad=False):
6061
"""Self-calibrate a dataset to an image.
@@ -64,9 +65,13 @@ def self_cal(obs, im, sites=[], method="both", pol='I', minimizer_method='BFGS',
6465
im (Image): the image to be calibrated to
6566
sites (list): list of sites to include in the self calibration.
6667
empty list calibrates all sites
67-
method (str): chooses what to calibrate, 'amp', 'phase', or 'both'
68-
minimizer_method (str): Method for scipy.optimize.minimize (e.g., 'CG', 'BFGS')
68+
6969
pol (str): which image polarization to self-calibrate visibilities to
70+
apply_singlepol (str): if calibrating to pol='RR' or pol='LL',
71+
apply solution only to the single polarization
72+
73+
method (str): chooses what to calibrate, 'amp', 'phase', or 'both'
74+
minimizer_method (str): Method for scipy.optimize.minimize (e.g., 'CG', 'BFGS')
7075
7176
pad_amp (float): adds fractional uncertainty to amplitude sigmas in quadrature
7277
gain_tol (float or list): gains that exceed this value will be disfavored by the prior
@@ -95,17 +100,21 @@ def self_cal(obs, im, sites=[], method="both", pol='I', minimizer_method='BFGS',
95100
"""
96101
if use_grad and (method=='phase' or method=='amp'):
97102
raise Exception("errfunc_grad in self_cal only works with method=='both'!")
98-
if pol not in ['I', 'Q', 'U', 'V', 'RR', 'LL']:
99-
raise Exception("Can only self-calibrate to I, Q, U, V, RR, or LL images!")
103+
100104
if pol in ['I', 'Q', 'U', 'V']:
101105
if obs.polrep != 'stokes':
102106
raise Exception("selfcal pol is a stokes parameter, but obs.polrep!='stokes'")
103107
im = im.switch_polrep('stokes', pol)
104-
elif pol in ['RR', 'LL', 'RRLL']:
108+
elif pol in ['RR', 'LL']:
105109
if obs.polrep != 'circ':
106110
raise Exception("selfcal pol is RR or LL, but obs.polrep!='circ'")
107111
im = im.switch_polrep('circ', pol)
108-
112+
else:
113+
raise Exception("Can only self-calibrate to I, Q, U, V, RR, or LL images!")
114+
115+
if apply_singlepol and obs.polrep!='circ':
116+
raise Exception("apply_singlepol must be False unless self-calibrating to 'RR' or 'LL'")
117+
109118
# V = model visibility, V' = measured visibility, G_i = site gain
110119
# G_i * conj(G_j) * V_ij = V'_ij
111120
if len(sites) == 0:
@@ -151,7 +160,7 @@ def self_cal(obs, im, sites=[], method="both", pol='I', minimizer_method='BFGS',
151160
if processes > 0: # run on multiple cores with multiprocessing
152161
scans_cal = np.array(pool.map(get_selfcal_scan_cal, [[i, len(scans), scans[i],
153162
im, V_scans[i], sites,
154-
obs.polrep, pol,
163+
obs.polrep, pol, apply_singlepol,
155164
method, minimizer_method,
156165
show_solution, pad_amp, gain_tol,
157166
debias, caltable, msgtype,
@@ -162,7 +171,7 @@ def self_cal(obs, im, sites=[], method="both", pol='I', minimizer_method='BFGS',
162171
for i in range(len(scans)):
163172
obsh.prog_msg(i, len(scans), msgtype=msgtype, nscan_last=i - 1)
164173
scans_cal[i] = self_cal_scan(scans[i], im, V_scan=V_scans[i], sites=sites,
165-
polrep=obs.polrep, pol=pol,
174+
polrep=obs.polrep, pol=pol, apply_singlepol=apply_singlepol,
166175
method=method, minimizer_method=minimizer_method,
167176
show_solution=show_solution,
168177
pad_amp=pad_amp, gain_tol=gain_tol,
@@ -208,7 +217,8 @@ def self_cal(obs, im, sites=[], method="both", pol='I', minimizer_method='BFGS',
208217
return out
209218

210219

211-
def self_cal_scan(scan, im, V_scan=[], sites=[], polrep='stokes', pol='I', method="both",
220+
def self_cal_scan(scan, im, V_scan=[], sites=[], polrep='stokes', pol='I', apply_singlepol=False,
221+
method="both",
212222
minimizer_method='BFGS', show_solution=False,
213223
pad_amp=0., gain_tol=.2, debias=True, caltable=False,
214224
use_grad=False):
@@ -223,6 +233,9 @@ def self_cal_scan(scan, im, V_scan=[], sites=[], polrep='stokes', pol='I', metho
223233
224234
polrep (str): 'stokes' or 'circ' to specify the polarization products in scan
225235
pol (str): which image polarization to self-calibrate visibilities to
236+
apply_singlepol (str): if calibrating to pol='RR' or pol='LL',
237+
apply solution only to the single polarization
238+
226239
method (str): chooses what to calibrate, 'amp', 'phase', or 'both'
227240
minimizer_method (str): Method for scipy.optimize.minimize
228241
(e.g., 'CG', 'BFGS', 'Nelder-Mead', etc.)
@@ -249,7 +262,7 @@ def self_cal_scan(scan, im, V_scan=[], sites=[], polrep='stokes', pol='I', metho
249262
sites = list(set(scan['t1']).union(set(scan['t2'])))
250263

251264
if len(V_scan) < 1:
252-
# This is not correct. Need to update to use polarization dictionary
265+
# TODO This is not correct. Need to update to use polarization dictionary
253266
uv = np.hstack((scan['u'].reshape(-1, 1), scan['v'].reshape(-1, 1)))
254267
A = obsh.ftmatrix(im.psize, im.xdim, im.ydim, uv, pulse=im.pulse)
255268
V_scan = np.dot(A, im.imvec)
@@ -343,10 +356,19 @@ def errfunc_grad(gpar):
343356
else:
344357
site_key = -1
345358

359+
# TODO: ANDREW - this has been changed
346360
# We will *always* set the R and L gain corrections to be equal in self calibration,
347361
# to avoid breaking polarization consistency relationships
348-
rscale = g_fit[site_key]**-1
349-
lscale = g_fit[site_key]**-1
362+
if apply_singlepol:
363+
if pol=='RR':
364+
rscale = g_fit[site_key]**-1
365+
lscale = np.ones(g_fit[site_key].shape)
366+
elif pol=='LL':
367+
rscale = np.ones(g_fit[site_key].shape)
368+
lscale = g_fit[site_key]**-1
369+
else:
370+
rscale = g_fit[site_key]**-1
371+
lscale = g_fit[site_key]**-1
350372

351373
# TODO: we may want to give two entries for the start/stop times
352374
# when a non-zero interval is used
@@ -357,23 +379,56 @@ def errfunc_grad(gpar):
357379
else:
358380
g1_fit = g_fit[g1_keys]
359381
g2_fit = g_fit[g2_keys]
360-
361382
gij_inv = (g1_fit * g2_fit.conj())**(-1)
362383

363384
if polrep == 'stokes':
385+
# gain factors
386+
g1_fit = g_fit[g1_keys]
387+
g2_fit = g_fit[g2_keys]
388+
gij_inv = (g1_fit * g2_fit.conj())**(-1)
389+
364390
# scale visibilities
365391
for vistype in ['vis', 'qvis', 'uvis', 'vvis']:
366392
scan[vistype] *= gij_inv
367393
# scale sigmas
368394
for sigtype in ['sigma', 'qsigma', 'usigma', 'vsigma']:
369395
scan[sigtype] *= np.abs(gij_inv)
396+
370397
elif polrep == 'circ':
371-
# scale visibilities
372-
for vistype in ['rrvis', 'llvis', 'rlvis', 'lrvis']:
373-
scan[vistype] *= gij_inv
374-
# scale sigmas
375-
for sigtype in ['rrsigma', 'llsigma', 'rlsigma', 'lrsigma']:
376-
scan[sigtype] *= np.abs(gij_inv)
398+
if apply_singlepol: #scale only solved polarization
399+
if pol=='RR':
400+
grr_inv = (g1_fit * g2_fit.conj())**(-1)
401+
gll_inv = np.ones(g1_fit.shape)
402+
grl_inv = (g1_fit)**(-1)
403+
glr_inv = (g2_fit.conj())**(-1)
404+
405+
elif pol=='LL':
406+
grr_inv = np.ones(g1_fit.shape)
407+
gll_inv = (g1_fit * g2_fit.conj())**(-1)
408+
grl_inv = (g2_fit.conj())**(-1)
409+
glr_inv = (g1_fit)**(-1)
410+
411+
# scale visibilities
412+
scan['rrvis'] *= grr_inv
413+
scan['llvis'] *= gll_inv
414+
scan['rlvis'] *= grl_inv
415+
scan['lrvis'] *= glr_inv
416+
417+
# scale sigmas
418+
scan['rrsigma'] *= np.abs(grr_inv)
419+
scan['llsigma'] *= np.abs(gll_inv)
420+
scan['rlsigma'] *= np.abs(grl_inv)
421+
scan['lrsigma'] *= np.abs(glr_inv)
422+
423+
else: #scale both polarizations
424+
gij_inv = (g1_fit * g2_fit.conj())**(-1)
425+
426+
# scale visibilities
427+
for vistype in ['rrvis', 'llvis', 'rlvis', 'lrvis']:
428+
scan[vistype] *= gij_inv
429+
# scale sigmas
430+
for sigtype in ['rrsigma', 'llsigma', 'rlsigma', 'lrsigma']:
431+
scan[sigtype] *= np.abs(gij_inv)
377432

378433
out = scan
379434

@@ -389,14 +444,14 @@ def get_selfcal_scan_cal(args):
389444
return get_selfcal_scan_cal2(*args)
390445

391446

392-
def get_selfcal_scan_cal2(i, n, scan, im, V_scan, sites, polrep, pol, method, minimizer_method,
447+
def get_selfcal_scan_cal2(i, n, scan, im, V_scan, sites, polrep, pol, apply_singlepol, method, minimizer_method,
393448
show_solution, pad_amp, gain_tol, debias, caltable, msgtype, use_grad):
394449
if n > 1:
395450
global counter
396451
counter.increment()
397452
obsh.prog_msg(counter.value(), counter.maxval, msgtype, counter.value() - 1)
398453

399-
return self_cal_scan(scan, im, V_scan=V_scan, sites=sites, polrep=polrep, pol=pol,
454+
return self_cal_scan(scan, im, V_scan=V_scan, sites=sites, polrep=polrep, pol=pol, apply_singlepol=apply_singlepol,
400455
method=method, minimizer_method=minimizer_method,
401456
show_solution=show_solution,
402457
pad_amp=pad_amp, gain_tol=gain_tol, debias=debias, caltable=caltable,

0 commit comments

Comments
 (0)