rcyndie / QuartiCal

CubiCal, but with greater power.
MIT License
0 stars 0 forks source link

Running a GKTB chain on LOFAR data #3

Open rcyndie opened 1 year ago

rcyndie commented 1 year ago

This is a MS (1 hr data, resolution: 4s, 488 channels) provided by Francesco. Please find below the phase and TEC plots he also provided (with the reference antenna 8).

Phase_plots TEC_plots

With a GKTB chain run on the data where, G - complex time variable K - delay time variable with an initial estimate (sol. time int. 10) T - TEC time variable with an initial estimate (sol. time int. 10) and B - complex frequency variable, we can get the following phase plots for the net effect of the chain image

There is "less smoothness" in the phase plots and as @JSKenyon pointed out, this might be because of the TEC and delay initial estimates are not aware of each other.

rcyndie commented 1 year ago

As per @JSKenyon's suggestion, a solution to this problem is running the GKTB chain in two parts, where there is a GK chain applied on the data and the resulting _partially correcteddata becomes the data column for the TB chain. Please find below the phase plots for the GK-net, TB-net and GKTB-net respectively.

image image image

(I redid the GKTB-net plots again as I had not referenced the antennas before calculating the net gains= GK-net*TB-net). The GKTB-net phase plot looks close to Francesco's. From this run of chain in parts, most of the effects are captured by the delay K and the remaining effects manifesting in the TB-net are mostly consistent with the ones in the B phase plots; please find below the phase plots for T and B respectively.

image image

The above results might be attributed to the narrow bandwidth of this observation, where there is an ambiguity between the delay and TEC effects.

rcyndie commented 1 year ago

The params TEC' from the solver tec_and_offset are calculated using the ionospheric phase delay definition $\phi = \exp (i2 \pi \text{TEC}')$, where $\text{TEC}' = \frac{-\kappa \text{TEC}}{c \nu}$. We can find the correct TEC using $\text{TEC} = \frac{- c \text{TEC}'}{10^{16} \kappa}$ in TECU. image While the shapes of the TEC graphs mostly agree in the above scale, the TEC values are smaller than Francesco's. This might be due to some rescaling constant.

landmanbester commented 1 year ago

This is excellent, well done. The shapes are almost identical so I guess its mostly the units. I suppose there may still be slight differences because it's not exactly the same solver and you are not solving for all the effects (at least I recall there being some additional ones in the paper Francesco shared with us). It would be interesting to look at the residual gains for each of the above. Let's discuss the way forward at the meeting tomorrow

rcyndie commented 1 year ago

Please also find below the phase plots from running a complex variable (say E) on solution time 1 and frequency 1 intervals. The variable E should essentially try its best to overfit the data and capture all (and any) possible effects with a higher d.o.f than the GKTB chain. Then the residual gains E - (GK*TB) should ideally contain mostly noise-like effects, confirming that the chain captured most of the effects. The phase plots are given for the E variable and the residual gains respectively.

image

image

The residual gains' plots are better demonstrated with a uniform colourmap (such as "viridis") rather than a cyclic one.

landmanbester commented 1 year ago

Interesting, there still seems to be a bit of structure in the residual gains. It's not obvious what it could be. Please will you prepare similar plots for GKB and GTB chains separately? That should give us an idea of how good the fit is for the different terms

rcyndie commented 1 year ago

The following is in support of the ambiguity between the delay and TEC effects mentioned earlier. If the ambiguity holds and we were to run a new chain GTKB in parts (GT and then KB), the results should tally with the GKTB chain.

The above results might be attributed to the narrow bandwidth of this observation, where there is an ambiguity between the delay and TEC effects.

Please find below the phase plots of the GT-net, KB-net and GTKB-net respectively. image image image

Although the phase plots match the GKTB chain, the TEC values are much higher than the ones obtained with the GKTB, as the clock (delay) corrections haven't been applied yet. image

E = complex variable with solution time 1 and frequency 1 intervals. Also consider the phase plots for the residual gains = E - (GT-net*KB-net),

image

JSKenyon commented 1 year ago

One other thing to note is that because of the contrived manner in which this chain is run, there may be some interplay between K and T which cannot be corrected.

Fixing this would require the initial estimate code to form corrected data on the fly. This may become somewhat nasty and some thought will need to be given to the memory footprint.

We can potentially do a little better by utilizing QuartiCal's ability to load and continue solving terms. To try this out, after you have solved for GK and TB as before, you can do one additional step where you solve for the entire GKTB chain. The difference is that you will need to disable the initial estimates and load the gain terms from disk. This is accomplished by setting e.g. G.load_from=path/to/previous/solution. Given that the K and T terms will already be close to solved, it should be safe to then iterate through the entire chain a couple of times. It would be interesting to see if it reduces the residual structure in the gains. Let me know if you need help with the additional configuration for this case @rcyndie!

rcyndie commented 1 year ago

The phase plots for the reference antenna are not strictly zero. Consider the amplitude and phase plots for the complex variable E and the GKB-net chain corresponding to the reference antenna (8).

image

The complex variable E is better able to fit "noise" with a higher degrees of freedom compared to the GKB chain. When set to the scale of the GKB chain, we can better visualise the "outliers" in the E plots. image

rcyndie commented 1 year ago

Please find the updated residual unwrapped phase plots:

  1. np.angle(E- GK*TB) image

  2. np.angle(E - GT*KB) image

  3. np.angle(GK*TB - GKB) image

4.np.angle(GK*TB - GTB) image

  1. np.angle(E - GKTB) (with loaded solutions and an extra solve) image

And the above can be viewed using a cyclic map

  1. np.angle(E- GK*TB) image

  2. np.angle(E - GT*KB) image

  3. np.angle(GK*TB - GKB) image

  4. np.angle(GK*TB - GTB) image

  5. np.angle(E - GKTB) (with loaded solutions and an extra solve) image

landmanbester commented 1 year ago

This is looking good @rcyndie. As discussed in the meeting, can you please make these plots by differencing the phases i.e.

np.angle(E) - np.angle(GKTB) 

and plot these with a cyclical colorbar. There are a number of additional experiments you can run:

1) GK (DATA, K initial estimate) followed by TB (CORRECTED_DATA, T initial estimate) then load GKTB and continue solve on DATA (no initial estimates) 2) GKTB (DATA, K initial estimate) and do 2 epochs 3) If 2 is successful repeat with B a constant in time complex 2x2 term.

Remember B = B(\nu) should always have an infinite time interval so you need to set the time chunking to 0. In addition to the gain residual plots you can have a look at making the chi2/dof plots with the surfchi2 app in surfvis

https://github.com/ratt-ru/surfvis

However, since the weights are not necessarily correct, we wouldn't expect a chi2/dof of one instead we need to compare these plots across experiments.

rcyndie commented 1 year ago

Following up on

https://github.com/rcyndie/QuartiCal/issues/3#issuecomment-1801468057,

please find below the difference in unwrapped phase plots calculated np.angle(a) - np.angle(b)

  1. np.angle(E) - np.angle(GK*TB) image

  2. np.angle(E) - np.angle(GT*KB) image

  3. np.angle(GK*TB) - np.angle(GKB) image

  4. np.angle(GK*TB) - np.angle(GTB) image

  5. np.angle(E) - np.angle(GKTB) (with loaded solutions and an extra solve) image

Also find the difference plots 11-15 when wrapped (and cyclic colourmaps).

  1. np.angle(E) - np.angle(GK*TB) image

  2. np.angle(E) - np.angle(GT*KB) image

  3. np.angle(GK*TB) - np.angle(GKB) image

  4. np.angle(GK*TB) - np.angle(GTB) image

  5. np.angle(E) - np.angle(GKTB) (with loaded solutions and an extra solve) image

rcyndie commented 1 year ago

Following up on the experiment ideas from

https://github.com/rcyndie/QuartiCal/issues/3#issuecomment-1801558645

For reference, also find the phase plots for a full GKTB chain without any initial estimates and 0a. solution time interval 10 image and 0b. solution time interval 1, image there is an increase in the flagged gain solutions.

Consider the difference in phase plots for solution time interval 10.

  1. np.angle(E - GKTB) (chain with only K init. est.) image

  2. np.angle(GKTB_cont - GKTB_kest) image

  3. np.angle(E) - np.angle(GKTB_kest) image

  4. np.angle(GKTB_cont) - np.angle(GKTB_kest) image

Consider the following residual phase plots when using a solution time interval of 1, as compared to solution time interval for plots 1 to 24,

  1. np.angle(E - GK*TB) (chain run in parts with initial estimates for K and T) image

  2. np.angle(E - GKTB) (chain run with loaded solutions from plot 25 and continued solve) image

  3. np.angle(E - GT*KB) (chain run in parts with initial estimates for T and K) << mostly sanity check and to be compared with plot 25 image

  4. np.angle(E - GKTB) (full chain with only K initial estimate) image

  5. np.angle(GKTB_cont - GKTB_kest) image

  6. np.angle(E) - np.angle(GK*TB) (chain run in parts with initial estimates for K and T) image

  7. np.angle(E) - np.angle(GKTB) (chain run with loaded solutions from plot 21 and continued solve) image

  8. np.angle(E) - np.angle(GT*KB) (chain run in parts with initial estimates for T and K) image

  9. np.angle(E) - np.angle(GKTB) (full chain with only K initial estimate) image

  10. np.angle(GKTB_cont) - np.angle(GKTB_kest) image

The GKTB run with loaded solutions and initial estimates from both K and T matches up closest with E, the overfitting complex variable in time and frequency. From plots 0a and 0b, it can be safely assumed that we cannot run a chain of the sort "GKTB" without expecting an increase in flagged gain solutions. For this data, if we were to rely on a difference of plots ph(a) - ph(b) (refer to plot 34), there is a "mildly" significant structure left when comparing a chain GKTB (loaded solutions + initial estimates on K and T) and a GKTB chain with initial estimate only on K, though the difference is more significant and dramatic in a ph(diff in gain solutions) (refer to plot 29).

rcyndie commented 1 year ago

Please find below the chi2 plots corresponding to the phases for the runs E, GKTB (chain with loaded solutions and init. est. on K and T), GKTB (chain with only K init. est.) and GKTB (chain with K init. est. and a complex 2x2 type for B but constant in time) using surfvis. I am sharing the plots corresponding to the diagonal correlations for the entire time and frequency chunks. NB: I might need to look at all the correlations for the GKTB chain run with K init. est. and 2x2 complex B. a: solution time interval 1 and b: solution time interval 10

35a. E (complex variable in time and frequency) image image

36a. GKTB (chain with loaded solutions and initial estimates with K and T) image image

36b. image image

37a. GKTB (chain with initial estimate only on K) image image

37b. image image

38a. GKTB (chain with initial estimate on K and a 2x2 complex B constant in time) image image

38b. image image

NB: The plots for the chains might not be the most accurate representation of the calibration runs, as these might be influenced by the supposedly "flagged" gain solutions still present in the data. Need to repeat the runs with solver_flags = TRUE and propagate_flags = TRUE. Coming shortly!

landmanbester commented 1 year ago

Ah, sorry. I should have sent you the command I meant for you to use. The above plots look like visibility waterfall plots made with the surfvis worker, I meant for you to make them with the surfchi2 worker. We can discuss it in the meeting shortly