Closed bennahugo closed 1 year ago
JSKenyon
commented
1 year ago Yep, this is a known problem which I have mentioned a few times already. It likely stems from the fact that the original interpolation code was proof of concept - it was never completely thought out. I am actively developing an improved version which properly handles phase (periodic quantities require special handling) and allows for correct interpolation of parameterised terms. Will let you know when it is complete!
JSKenyon
commented
1 year ago It would be worth revisiting these experiments with the v0.2.0-dev branch (no urgency though). The interpolation has been overhauled, and reference antennas are honoured on diagonal terms. If the problem persists, I will go on another bug hunt.
bennahugo
commented
1 year ago I will try to as soon as I have time again. It is worth doing a direct transfer comparison with casa over a wide swath to see how the scatter on the corrected data compare on a set of calibrator fields.
On Thu, May 25, 2023 at 11:20 AM JSKenyon @.***> wrote:
It would be worth revisiting these with the v0.2.0-dev branch now. The interpolation has been overhauled, and reference antennas are honoured on diagonal terms. If the problem persists, I will go on another bug hunt.
— Reply to this email directly, view it on GitHub https://github.com/ratt-ru/QuartiCal/issues/251#issuecomment-1562570708, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB4RE6SY3LDIRGEUXKJUQHDXH4P5BANCNFSM6AAAAAAXU6W7GM . You are receiving this because you authored the thread.Message ID: @.***>
Benjamin Hugo
PhD. student, Centre for Radio Astronomy Techniques and Technologies Department of Physics and Electronics Rhodes University
Junior software developer Radio Astronomy Research Group South African Radio Astronomy Observatory Black River Business Park Observatory Cape Town
bennahugo
commented
1 year ago @JSKenyon Retrying this with the latest branch on 1934 -> 3C286 transfer on S0 data with #fe1a23d96828b. I restore flags with casa flagmanager prior to switching software.
CASA <17>: gaincal(vis='s0_offsets.flagged.avg.uncorrected.ms', caltable='1934.K', refant='2', solint='64s'
...: , field='1', gaintype='K')
CASA <19>: gaincal(vis='s0_offsets.flagged.avg.uncorrected.ms', caltable='1934.Gp', refant='2', solint='64s
...: ', field='1', gaintype='G', calmode='p', gaintable='1934.K')
CASA <20>: gaincal(vis='s0_offsets.flagged.avg.uncorrected.ms', caltable='1934.Ga', refant='2', solint='64s
...: ', field='1', gaintype='G', calmode='a', gaintable=['1934.K','1934.Gp'])
CASA <23>: bandpass(vis='s0_offsets.flagged.avg.uncorrected.ms', caltable='1934.B', field='1', refant='2',
...: solint='inf', combine='scan', gaintable=['1934.K', '1934.Gp', '1934.Ga'])
CASA <24>: applycal(vis='s0_offsets.flagged.avg.uncorrected.ms', field='0,1', gaintable=['1934.K', '1934.Gp
...: ', '1934.Ga', '1934.B'])yields primary corrected data:
and 3C286 corrected data:
As far as I can tell the following should be equivalent to solving per scan and then across scan boundary as is done in casa
### KG per scan
goquartical \
input_ms.path=s0_offsets.flagged.avg.uncorrected.ms \
input_ms.data_column=DATA \
input_ms.time_chunk=100000000s \
input_ms.select_fields="[1]" \
input_ms.group_by="['SCAN_NUMBER',FIELD_ID,DATA_DESC_ID]" \
input_model.recipe=MODEL_DATA \
solver.terms="['K','Gp','Ga']" \
solver.iter_recipe="[0,20,20]" \
solver.reference_antenna="2" \
output.products="[]" \
output.columns="[]" \
output.overwrite=False \
Gp.type=phase \
Gp.time_interval=64s \
Gp.freq_interval=0 \
Ga.type=amplitude \
Ga.time_interval=64s \
Ga.freq_interval=0 \
K.type=delay \
K.time_interval=64s \
K.freq_interval=0 \
K.initial_estimate=True
### B inf, combine="scan"
goquartical \
input_ms.path=s0_offsets.flagged.avg.uncorrected.ms \
input_ms.data_column=DATA \
input_ms.time_chunk=100000000s \
input_ms.select_fields="[1]" \
input_ms.group_by="[FIELD_ID,DATA_DESC_ID]" \
input_model.recipe=MODEL_DATA \
solver.terms="['K','Gp','Ga','B']" \
solver.iter_recipe="[0,0,0,150]" \
solver.reference_antenna="2" \
output.products="[]" \
output.columns="[]" \
output.overwrite=False \
output.gain_directory=gains.qc.B \
Gp.type=phase \
Gp.time_interval=64s \
Gp.freq_interval=0 \
Gp.load_from=gains.qc/Gp \
Ga.type=amplitude \
Ga.time_interval=64s \
Ga.freq_interval=0 \
Ga.load_from=gains.qc/Ga \
B.type=diag_complex \
B.time_interval=1000000000s \
B.freq_interval=1 \
K.type=delay \
K.time_interval=64s \
K.freq_interval=0 \
K.initial_estimate=False \
K.load_from=gains.qc/K
### applycal KGB to both fields 0,1
goquartical \
input_ms.path=s0_offsets.flagged.avg.uncorrected.ms \
input_ms.data_column=DATA \
input_ms.time_chunk=100000000s \
input_ms.select_fields="[0,1]" \
input_ms.group_by="[FIELD_ID,DATA_DESC_ID]" \
input_model.recipe=MODEL_DATA \
solver.terms="['K','Gp','Ga','B']" \
solver.iter_recipe="[0,0,0,0]" \
solver.reference_antenna="2" \
output.products="[corrected_data]" \
output.columns="['CORRECTED_DATA']" \
output.overwrite=False \
output.gain_directory=gains.qc.sec.transferonly \
Gp.type=phase \
Gp.time_interval=64s \
Gp.freq_interval=0 \
Gp.load_from=gains.qc/Gp \
Ga.type=amplitude \
Ga.time_interval=64s \
Ga.freq_interval=0 \
Ga.load_from=gains.qc/Ga \
B.type=diag_complex \
B.time_interval=1000000000s \
B.freq_interval=1 \
B.load_from=gains.qc.B/B \
K.type=delay \
K.time_interval=64s \
K.freq_interval=0 \
K.initial_estimate=False \
K.load_from=gains.qc/KThis yields a much wider spread on the transferred source:
Not sure where the issue is here, but it looks like the transfer calibration is not quite working yet, unless I made a blunder
bennahugo
commented
1 year ago Looks like delays going absolutely bos
bennahugo
commented
1 year ago I should add that the corrected phase on the 1934 where the delay is solved looks almost as good as CASA so it is probably a failure in interpolation?
bennahugo
commented
1 year ago @JSKenyon as discussed the interpolation onto 0,1 simultaneous and onto 0 and 1 as separate steps appears to be about equal
Here is the primary for reference - the delays look reasonable
3C286 is as above.
bennahugo
commented
1 year ago Alright @JSKenyon and I tracked this down. Note this issue @landmanbester. This is due to solutions being linearly extrapolated in time. The observation is set up as follows:
2023-07-18 07:49:56 INFO listobs::ms::summary+ 31-Jan-2023/00:37:58.5 - 00:43:02.5 1 0 J1331+3030 1584 [0] [101] [UNKNOWN]
2023-07-18 07:49:56 INFO listobs::ms::summary+ 00:59:34.5 - 01:04:22.5 12 0 J1331+3030 1584 [0] [96] [UNKNOWN]
2023-07-18 07:49:56 INFO listobs::ms::summary+ 01:20:54.6 - 01:25:50.6 23 0 J1331+3030 1584 [0] [98.7] [UNKNOWN]
2023-07-18 07:49:56 INFO listobs::ms::summary+ 01:43:26.6 - 01:46:14.6 34 1 J1939-6342 1056 [0] [84] [CALIBRATE_BANDPASS,CALIBRATE_FLUX]
2023-07-18 07:49:56 INFO listobs::ms::summary+ 01:47:42.7 - 01:52:30.7 35 0 J1331+3030 1584 [0] [96] [UNKNOWN]
2023-07-18 07:49:56 INFO listobs::ms::summary+ 02:10:14.7 - 02:13:02.7 46 1 J1939-6342 1056 [0] [84] [CALIBRATE_BANDPASS,CALIBRATE_FLUX]
2023-07-18 07:49:56 INFO listobs::ms::summary+ 02:14:30.7 - 02:19:26.7 47 0 J1331+3030 1584 [0] [96.6] [UNKNOWN]
2023-07-18 07:49:56 INFO listobs::ms::summary+ 02:37:02.8 - 02:39:58.8 58 1 J1939-6342 1056 [0] [87.2] [CALIBRATE_BANDPASS,CALIBRATE_FLUX]
2023-07-18 07:49:56 INFO listobs::ms::summary+ 02:41:18.8 - 02:46:14.8 59 0 J1331+3030 1584 [0] [98.7] [UNKNOWN]
2023-07-18 07:49:56 INFO listobs::ms::summary+ 03:03:58.8 - 03:06:54.9 70 1 J1939-6342 1056 [0] [88] [CALIBRATE_BANDPASS,CALIBRATE_FLUX]
2023-07-18 07:49:56 INFO listobs::ms::summary+ 03:08:22.9 - 03:13:18.9 71 0 J1331+3030 1584 [0] [98.7] [UNKNOWN]
2023-07-18 07:49:56 INFO listobs::ms::summary+ 03:31:02.9 - 03:33:58.9 82 1 J1939-6342 1056 [0] [88] [CALIBRATE_BANDPASS,CALIBRATE_FLUX]
2023-07-18 07:49:56 INFO listobs::ms::summary+ 03:35:34.9 - 03:40:22.9 83 0 J1331+3030 1584 [0] [96] [UNKNOWN]
2023-07-18 07:49:56 INFO listobs::ms::summary+ 03:58:15.0 - 04:01:11.0 94 1 J1939-6342 1056 [0] [87.8] [CALIBRATE_BANDPASS,CALIBRATE_FLUX]
2023-07-18 07:49:56 INFO listobs::ms::summary+ 04:02:47.0 - 04:07:35.0 95 0 J1331+3030 1584 [0] [96] [UNKNOWN]
2023-07-18 07:49:56 INFO listobs::ms::summary+ 04:25:27.0 - 04:28:23.1 106 1 J1939-6342 1056 [0] [88] [CALIBRATE_BANDPASS,CALIBRATE_FLUX]
2023-07-18 07:49:56 INFO listobs::ms::summary+ 04:29:59.1 - 04:34:55.1 107 0 J1331+3030 1584 [0] [98.7] [UNKNOWN]
2023-07-18 07:49:56 INFO listobs::ms::summary+ 04:52:47.1 - 04:55:43.1 118 1 J1939-6342 1056 [0] [88] [CALIBRATE_BANDPASS,CALIBRATE_FLUX]
2023-07-18 07:49:56 INFO listobs::ms::summary+ 04:57:27.1 - 05:02:15.1 119 0 J1331+3030 1584 [0] [96] [UNKNOWN]
2023-07-18 07:49:56 INFO listobs::ms::summary+ 05:20:23.2 - 05:23:11.2 130 1 J1939-6342 1056 [0] [84] [CALIBRATE_BANDPASS,CALIBRATE_FLUX]
2023-07-18 07:49:56 INFO listobs::ms::summary+ 05:25:03.2 - 05:29:51.2 131 0 J1331+3030 1584 [0] [96] [UNKNOWN]
Ie. 3C286 is observed prior to the first delay solution on 1934. When extrapolated over a wide span the solutions, in K, are completely invalid (you can imagine transferring the delays across many hours - they are markedly stable along with the passband. I believe CASA will use nearest neighbour for the edges of the domain so we don't see this issue.
When plotting the Quartical corrected data by scan you can clearly see the effect of extrapolating with a gradient show up:
In CASA these are, for reference
Note scan 35 which falls inside the time domain of the interpolation. This is almost 1:1 to my eye, so I'm happy the interpolation itself works inside the domain - only the edgecases (out of domain / flagged solutions) need to be handled @JSKenyon
bennahugo
commented
1 year ago I think the way to do this is perhaps to use https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html#scipy-interpolate-interp1d with bounds_error set and in a try catch pick the nearest unflagged value?
JSKenyon
commented
1 year ago With the changes in #285, QC produces the following plots for the same scans as posted by @bennahugo above. These are very consistent (by eye) with the casa results. It is also worth noting that it may be possible to get slightly better results with different QC settings - these results were generated with a similar procedure to casa but don't use the full extent of QC functionality.
JSKenyon
commented
1 year ago As a follow up, I tried doing something similar in a smaller number of steps. These configs (configs.tar.gz) produce the following results which may be fractionally better:
bennahugo
commented
1 year ago I believe this is now working as expected.
I suggest we use that dataset to build a transfer calibration acceptance test.
--- Here is a bit of a story for comparison sake - sorry for the read ---
I'm trying to execute a recipe analogous to a casa-style transfer calibration without achieving good results from interpolation. Something is markedly off in phase transfer.
The goal is to get absolute leakage corrections for the reference antenna I picked in an earlier reduction in CASA to get error bar estimates for linear angle offsets ala Hales 2017 of the real part of the leakage solution on the reference antenna that was artificially set to 0 by casa (I suspect this is minimal -- things are solidly bound by a 1rad/m2 error from the best ionospheric corrections in the business...)
I need scale and bandpass transferred from 1934 because we can only calibrate bandpasses with unpolarized calibrators (we will never know the ionosphere and intrinsic angles well enough to do this with polarized sources without severe biases in the case of linear feed bases)
For the moment I forgo delays and just mop them into the bandpass in the transfer not to overly complicate interpolation. They should be relatively stable as they are cable delays - modulo geometric errors due to residual antenna position errors on MK.
I fill the model column (with setjy to Stevens Reynolds 2016 for PKS 1934 - the primary) and fill the model column for 3C286 to unity using CASA clearcal.
(This is S0 band data on MK with only calibrators - I'm calibrating calibrators. I solve across scan boundary to give me 'inf' combine='scan' capability. The scans are closeby temporally)
Primary looks good
Then we transfer the phase and everything over to inspect
Transfer phases are completely non-sensical. Scale is very roughly where I suspect it to be though
The equivalent GpGaB transfer-only operation in CASA produces something much more sensible (phases are not perfect ofc due to the separation distance, but things are not nearly as bad, even with delays completely foregone)
Thus: quartical interpolation seems not to be working as expected
This is not a huge issue for me per say due to my observation not being imaging, but I suspect this will be a problem for the final calibration step onto a target for a standard imaging observation. This fixes the phase very roughly (the bandpass mopped delays causes a bit of a classic residual delay spread)
If I include the delay (by lagspace method on the primary) during the solve for the bandpass but discard it and resolve for it towards the secondary (ie apply GaB and solve KGp) then I get something that does not look much better. The equivalent operation in CASA gives something with equal spread across the entire band
If I make K solvable to 250 iterations things improve markedly :)
This is good news and bad news - the phase transfer obviously needs a fix