Notebooks

[NicoRenaud]

Dear Hera Team, thanks for all your work and for making all your tools available. I'm trying to execute the omnical_convergence.ipynb notebooks that is in the scripts/notebooks folder but it seems to be out of sync with the library. Most of the issues seem to be related with the om.RedSol class that I think was introduced after the last update of the notebook.

I've managed to fix a few of these issues but I am not sure if what I'm doing is correct.
Would it be possible for you to update the notebooks ?

Thanks again

[jsdillon]

Hi Nico,

As you may have noticed, those notebooks haven't been updated in 5 years but the interfaces have changed considerably in that time to adapt to the needs of HERA analysis. In particular, the wrapper around the C of the old omnical (written by Jeff Zheng) has been completely deprecated. So in short, I don't think the notebook/memo can be trivially updated.

If you really want to run this memo, I'd suggest checking out an old commit of this repo and trying that. I'm not sure how old you have to go and what version of related software you'd also need to check out to make that work. Might be very annoying to actually get it to run without errors, so fair warning. In general, this repo is meant more as a way for HERA team members to collaborate and a version-controlled record of how HERA analysis was performed, rather than a public resource for the broader astronomical community (which comes with certain expectations about interface stability, deprecation warnings, release notes, etc.) We just don't have the person-power to do that right now. That said, if our codebase can be helpful to others, we certainly want to share! Perhaps you can tell me more about what you're trying to do and I can point you in the right direction?

One thing I'll also note is that I think that that's notebook description of omnical as "linsolve based" is a bit misleading. It does use that code, but doesn't actually use any of the solvers in that code---just the infrastructure for setting up systems of equations. You can read more about what omnical is doing here: https://arxiv.org/abs/2003.08399

[NicoRenaud]

Hi Josh,

thanks a lot for your quick and detailed answer ! I fully understand the issue of not having enough person-power. I was expecting this to be the issue but I thought I would simply ask what was the status of the notebooks.

I'm working with people from ASTRON in the Netherlands, and in particular Chris Broekema @broekema. We are exploring the possibility of integrating quantum computing in various radio astronomy pipelines. We are not expecting any advantages nor speed up given the current status of the technology, but we want to illustrate the role that quantum computers could play in the future.

One of the option we have explored is to simply replace the linear solvers used in redcal by quantum linear solvers. We have recently finished implementing one quantum solver that can be deployed on near term quantum chips (IBM-Quantum machines) and we are now starting to integrate it in hera_cal via linsolve

We are finishing the paper and I would like to have an integration use case where we use our quantum solvers on real data. I was hoping to use the examples in scripts/notebooks/ as a starting point for that.

So if you have an example where hera_cal is used to perform calibration on a small amount of data and if you are willing to point me in its direction that would be very much appreciated !

Thanks a lot !

[jsdillon]

That's super cool! Sounds very futuristic, but I'm glad someone is exploring the idea of using quantum computers for radio astronomy.

So the big question here is precisely which algorithms you want to compare. There's:

• logcal: a non-iterative but biased solver for the system of equations
• lincal: an iterative, unbiased solver, but much slower approach that uses a Gauss-Newton algorithm. It's O(N^3) where N is the number of antennas.
• omnical: an alternative iterative algorithm that takes many more steps (each one much, much faster) to achieving what lincal does in O(N^2) time using damped fixed-point iteration (check out Figure 2 of https://arxiv.org/pdf/2003.08399.pdf)

In our current calibration of HERA data, we're only using omnical, which barely uses any of the guts of linsolve. The actual damped fixed-point solver is here:

hera_cal/hera_cal/redcal.py

Line 847 in d0e61fa

def solve_iteratively(self, conv_crit=1e-10, maxiter=50, check_every=4, check_after=1,

[NicoRenaud]

Hi Josh,

I was targeting the firstcal stage. But I think 'logcal' and 'lincal' could also be good options to use the quantum solver. I haven't tried that yet though. So if there is somewhere a small example/notebook that uses these solvers on a small dataset that would be awesome !

Thanks for your time !

[jsdillon]

Hi Nico,

Let me see if I can throw together a script/notebook that does what you want using more up-to-date code.

-Josh

[jsdillon]

Hi Nico,

Here's a basic script that demonstrates redundant-baseline calibration. It generates a 61-element array with reasonable true gains and foregrounds, then performs calibration and shows that you can get out what you put in (after fixing some of the degeneracies of redundant-baseline calibration).

import numpy as np
from hera_sim.antpos import hex_array, linear_array
from hera_sim.vis import sim_red_data
from hera_cal import redcal, abscal, apply_cal
from hera_cal.utils import split_pol
from hera_sim.sigchain import Bandpass
from hera_sim.foregrounds import DiffuseForeground, PointSourceForeground

# Generate array layout
antpos = hex_array(5, split_core=False, outriggers=0)
reds = redcal.get_reds(antpos, pols=['ee'], pol_mode='1pol')

# Define frequencies 
freqs = np.linspace(120e6, 180e6, 100) # in Hz
df = np.median(np.diff(freqs))

# Gains with a bandpass model and cable delays bandpass with cable delays
bp = Bandpass(gain_spread=.1, dly_rng=(-20, 20))
gains = {(ant, 'Jee'): gain[None, :] for ant, gain in bp(freqs / 1e9, list(antpos.keys())).items()}

vis = {red[0]: PointSourceForeground()(np.array([0]), freqs / 1e9,  antpos[red[0][1]] - antpos[red[0][0]]) for red in reds}

# Build RedSol object with true solutions
sol_true = redcal.RedSol(reds, gains=gains, vis=vis)

# Create uncalibrated data
data = {bl: sol_true[red[0]] * sol_true[bl[0], 'Jee'] * sol_true[bl[1], 'Jee'].conj() for red in reds for bl in red}

# Perform redundant-baseline calibration
rc = redcal.RedundantCalibrator(reds)
meta_fc, sol_fc = rc.firstcal(data, freqs)
meta_logcal, sol_logcal = rc.logcal(data, sol0=sol_fc)
meta_omnical, sol_omnical = rc.omnical(data, sol_logcal, maxiter=500, check_after=50, check_every=10)

# Assert that chi^2 is 0 to high precision
np.testing.assert_almost_equal(meta_omnical['chisq'], 0, decimal=8)

# Perform absolute calibration to get close enough that remove_degen() can fix degneracies of redundant-baseline calibration
dly_slope_gains = abscal.delay_slope_lincal(sol_true.vis, sol_omnical.vis, antpos, df=df, f0=freqs[0], medfilt=False,
                                            assume_2D=True, verbose=False, return_gains=True, gain_ants=sol_omnical.gains.keys())
sol_omnical.gains = {antpol : g * dly_slope_gains[antpol] for antpol, g in sol_omnical.gains.items()}
apply_cal.calibrate_in_place(sol_omnical.vis, dly_slope_gains)

phase_slope_gains = abscal.global_phase_slope_logcal(sol_true.vis, sol_omnical.vis, antpos, reds=reds, solver='linfit', assume_2D=True,
                                               verbose=False, return_gains=True, gain_ants=sol_omnical.gains.keys())
sol_omnical.gains = {antpol : g * phase_slope_gains[antpol] for antpol, g in sol_omnical.gains.items()}
apply_cal.calibrate_in_place(sol_omnical.vis, phase_slope_gains)    

abscal_meta, delta_gains = abscal.complex_phase_abscal(sol_omnical.vis, sol_true.vis, sol_omnical.reds, 
                                                       list(sol_omnical.vis.keys()), list(sol_true.vis.keys()))
sol_omnical.gains = {antpol : g * delta_gains.get(antpol, 1) for antpol, g in sol_omnical.gains.items()}
apply_cal.calibrate_in_place(sol_omnical.vis, delta_gains)     
sol_fc.remove_degen(degen_sol=sol_true)

# Finally fix degenerate modes to exactly match true values
sol_omnical.remove_degen(degen_sol=sol_true)

# Assert visibility solutions and gains are correct to high precision
for key in sol_omnical:
    np.testing.assert_almost_equal(sol_omnical[key], sol_true[key], decimal=8)

Right now, it's noise-free, so another thing you could try is adding noise to the data. You won't be able to show that you get the exact right answer, but you could show that the normalized chi^2 (see https://github.com/HERA-Team/hera_cal/blob/d0e61fa6d1c8aad5f4e9cae2bdbb99ba83830fca/hera_cal/redcal.py#L1479C1-L1479C1) has the right distribution. There's more about that in the paper I mentioned above.

Let me know if this is the kind of thing you were looking for.

[NicoRenaud]

Hey Josh,

thanks a lot for your time and efforts ! That's amazing !
I'll dive into the scripts and see how I can use it for our purpose. It works already great and I can plug in our quantum redundant solver so it's already a great start.

Thanks again for all your help