Common beam methods

[e-koch]

The continuation of #62.

I've filled out the common beam methods to a point where they can be used. There remains some small quirks that I don't understand, but they result in almost negligible changes to the common beam.

Methods

1. Pairwise ellipse comparison - common_2beams. This is based on the CASA implementation and has an exact analytical result only for sets of 2 beams. Included in the tests is the same testing cases in the CASA tests, though I found a few cases where the target common beams were not the optimal solution. These are highlighted in the test cases.

2. @low-sky's optimization method - common_manybeams_opt. Write-up here with an implementation here. I made some small changes to the objective function, but there are many cases where the opt.minimize can't find the global minimum. A brute-force search with opt.brute works better, but it isn't obvious to me how to define a small enough parameter space to make this a quick operation. For now, common_manybeams_opt raises NotImplementedError.

3. Minimum volume ellipse (MVE) from the Khachiyan Algorithm - common_manybeams_mve. I adapted the implentation from here, which comes from the original here. This method is close to the optimization method, but it requires a set of bounding points to fit the MVE to. To do this, a number of sample points are created for each ellipse (default to 200 each). The bounding points are then identified with a convex hull, which requires scipy as a dependency. The minimum enclosing ellipse is fit to the edge points of the convex hull. The Khachiyan Algorithm find the optimal solution by maximizing det(A), where A is the center-form of the ellipse. The solution converges to a given error tolerance (default 1e-4). Decreasing the tolerance by 1/N requires N times more iterations to converge, and after various tests, a tolerance of ~1e-4 seems to be a good comprise between time and precision.

I've found that this method slightly underestimates the MVE, which leads to a common beam that marginally cannot be deconvolved by all beams in the set. There is inherent uncertainty due to the finite sampling of the ellipse edges and the error tolerance, but the underestimate remains when sampling with more edge points or requiring a lower tolerance. I am not sure where this arises from, but the difference between the MVE and the actual common beam is small. I introduced a fudge factor epsilon that increases the radii of all edge samples by (1 + epsilon). Setting epsilon=5e-4 gives a common beam that is deconvolvable for all of the test cases and within ~0.3% of the true value. Note that the CASA implementation of method 1 also includes a small fudge factor in cases where the resulting common beam is not deconvolvable.

Beams.common_beam() uses method 1 for any set of two ellipses and method 3 for all other cases. The tests include consistency checks between the two methods.

Other changes

• Added Beam.iscircular(rtol=1e-6) to check for circular beams to within some tolerance range. The PA angle checks in Beam.__eq__ are now only performed if the beam is non-circular. The tests were also updated for checks with circular beams.

[e-koch]

@low-sky - You're not coming up on the request review list, but you should also review this PR, please.

[keflavich]

@e-koch Keep good notes on the process of developing this algorithm. It will go into the radio-astro-tools paper when we write it up (...later this year?)

[low-sky]

FWIW, the optimization method really does seem to require the convex optimization for a stable solution. It's a classic problem in the convex optimization world, so I think it's reasonable that this wall would be encountered. I'm not really an expert here but it may be possible to do a stripped-down implementation of the single problem we have to solve.