http://docs.julialang.org/en/latest/stdlib/math.html#Base.cis

Didn't know this function existed tbh.

At the moment TTCal prints a welcome banner every time it's run, but doesn't actually describe what it is currently doing at all.

Currently the calibration routine unnecessarily uses too much memory. A full band LWA measurement set just barely squeeks by with 93% memory use on my machine.

Counting in units of the size of the data:

• 1 for the data itself
• 1 for the model
• 2 for the data packed into square matrices
• 2 for the model packed into square matrices

For a 5 GB dataset we're looking at 30 GB of memory which is way too much.

The last two points are where we can improve because those matrices can be computed as we iterate through frequency channels instead of doing it all at once at the beginning. Note also that the factor of 2 comes because we populate both the upper and lower triangular parts.

I don't think this problem affects the master branch.

Currently Travis is consistently timing out on the calibration tests because they go >10 min without outputting anything. Notably it is not hanging. If I insert some print statements you can see that it's still running, just incredibly slowly. I feel like I should have noticed this locally though.

I've managed to rule out the latest changes to CasaCore.jl. So there's either something sinister in this commit or something about Julia v0.5 is slowing things down. Inlining changes maybe?

This issue seems relevant: JuliaLang/julia#17476

related to #28

bad antennas and frequency channels should be detected during calibration and flagged

I have a prototype using NLopt to do the position fitting. It is so much more robust than the janky piece of code I wrote myself. I'm currently using the COBYLA algorithm for its support of equality constraints.

    count = 0
    start = [0.0, 0.0, 0.0]
    objective(x, g) = (count += 1; fitvis_nlopt_objective(x, g, data, flags, meta, direction, uvw))
    constraint(x, g) = fitvis_nlopt_constraint(x, g, direction)

    opt = Opt(:LN_COBYLA, 3)
    equality_constraint!(opt, constraint)
    max_objective!(opt, objective)
    initial_step!(opt, sind(10/60))
    xtol_abs!(opt, tolerance)
    flux, offset, _ = optimize(opt, start)
    @show count flux offset _

The equality constraint is crucial for optimizations on the unit sphere. This is because I would rather solve for the components of the unit vector than the spherical coordinates (theta, phi) to avoid the coordinate discontinuity. However you need to enforce the constraint that the magnitude of the unit vector is one.

NLopt does have an augmented Lagrangian method that helps you attach equality constraints to methods that don't natively support them. I got this to work with the Nelder-Mead method as an example but I couldn't figure out the right incantation that would get a reasonable stopping criterion (it was taking ~10,000 iterations to converge irrespective of the tolerance).

    # NB I wasn't really able to figure out what needs to be attached to `opt`
    # versus what needs to be attached to `local_opt`.
    opt = Opt(:AUGLAG, 3)
    local_opt = Opt(:LN_NELDERMEAD, 3)
    local_optimizer!(opt, local_opt)
    equality_constraint!(opt, constraint)
    max_objective!(opt, objective)
    initial_step!(local_opt, sind(10/60))
    xtol_abs!(opt, sqrt(tolerance))
    flux, offset, _ = optimize(opt, start)
    @show count flux offset _

Edit: I also noticed that if I set a maximum number of objective evaluations, and the optimization was terminated for exceeding this maximum, then NLopt would give me the original starting point instead of its current best guess.

I just noticed that ttcal-0.2.0 errors out if you attempt to shave a near-field source using a source model flux where the minimum specified frequency does not go below the lowest frequency of the measurement set. Using a sources.json file that looks like your example in the TTCal docs,
[
{
"name": "Noisy Power Lines",
"sys": "WGS84",
"long": -118.31478,
"lat": 37.14540,
"el": 1226.709,
"rfi-frequencies": [3.0e7, 3.75e7, 4.5e7, 5.25e7, 6.0e7],
"rfi-I": [1.0, 2.0, 3.0, 4.0, 5.0]
}
]

will produce the following error when used to shave subband 0:

$ ttcal-0.2.0 shave --input 00-T1al.ms --sources sources.json --minuvw 10 --maxiter 30 --tolerance 1e-4
Running shave on 00-T1al.ms
ERROR: LoadError: BoundsError: attempt to access 5-element Array{TTCal.StokesVector,1}:
TTCal.StokesVector(1.0,0.0,0.0,0.0)
TTCal.StokesVector(2.0,0.0,0.0,0.0)
TTCal.StokesVector(3.0,0.0,0.0,0.0)
TTCal.StokesVector(4.0,0.0,0.0,0.0)
TTCal.StokesVector(5.0,0.0,0.0,0.0)
at index [0]
in genvis_onesource! at /opt/astro/mwe/ttcal-0.2.0/julia-packages/v0.4/TTCal/src/genvis.jl:142
in genvis at /opt/astro/mwe/ttcal-0.2.0/julia-packages/v0.4/TTCal/src/genvis.jl:25
in shave! at /opt/astro/mwe/ttcal-0.2.0/julia-packages/v0.4/TTCal/src/peel.jl:86
in run_shave at /opt/astro/mwe/ttcal-0.2.0/julia-packages/v0.4/TTCal/src/commandline.jl:267
in main at /opt/astro/mwe/ttcal-0.2.0/julia-packages/v0.4/TTCal/src/commandline.jl:199
in include at ./boot.jl:261
in include_from_node1 at ./loading.jl:320
in process_options at ./client.jl:280
in _start at ./client.jl:378
while loading /opt/astro/mwe/ttcal-0.2.0/julia-packages/v0.4/TTCal/src/ttcal.jl, in expression starting on line 19

If I change the first element of “rfi-frequencies” from 3.0e7 —> 2.0e7, it runs without a hitch. Also, this error does not happen for measurement sets where the top of the subband is higher than the highest frequency specified in the source model — i.e. running shave on subband 21 with the same sources.json above does not produce an error.

This isn’t impacting me (I can easily change my sources.json file!), but just wanted to give a heads up!

• Marin

I should be able to eek out some more performance by using StructsOfArrays.jl to SIMD vectorize the inner calibration loop.

I need to double check the inner calibration step for polcal. Specifically, it computes numerator/denominator. Is this the correct operation? If so, I need to define /(::JonesMatrix,::JonesMatrix). Otherwise I should change it to denominator\numerator.

• test iterative fringe calculation against the naive method
• test flagdata
• test bandpass with no noise
• test bandpass with noise added to the autocorrelations
• test bandpass with noise added to everything
• test polcal
• test applycal
• test source fitting
• test selfcal
• test the TTCal.run(...) interface in all cases!
• test the command line interface in all cases

Add a more descriptive error message for people attempting to use an out of date sources.json.

The error comes from readsources looking for the Stokes parameters as keys.

At the moment I set the gains to unity and flag the channel. This is somewhat annoying in practice.

Edit: especially annoying if the gains were somewhat close to converging

Peeling is currently very sensitive to unmodeled flux. Especially RFI sitting on the horizon.

@marinanderson You should probably double check my math, but I think I found an error in the stefcal step for polarized calibration here:

I think the update step actually needs to be conj(denominator\numerator) instead of (denominator\numerator)' (which does a conjugate transpose).

needs to be (plasma frequency / frequency)^2

edit: the implementation is correct though