predict_checks fails for linear parallel hand only data

[bennahugo]

With fixes applied for #31 (spi set to [0,0,0,0] and reffreq set to 0.0) for flat spectrum prediction. I get

2020-09-26 10:40:03 | INFO | predict:parse_sky_models | Source groups/clusters for /home/hugo/workspace/DEEP2_tagged.lsm.html:
  DIE     : 0 point source/s, 1364 Gaussian source/s
2020-09-26 10:40:03 | ERROR | goquartical:<module> | An error has been caught in function '<module>', process 'MainProcess' (28998), thread 'MainThread' (140358578743104):
Traceback (most recent call last):

  File "/usr/lib/python3.6/runpy.py", line 193, in _run_module_as_main
    "__main__", mod_spec)
                └ ModuleSpec(name='__main__', loader=<_frozen_importlib_external.SourceFileLoader object at 0x7fa7c72464e0>, origin='/home/hugo...

  File "/usr/lib/python3.6/runpy.py", line 85, in _run_code
    exec(code, run_globals)
         │     └ {'__name__': '__main__', '__doc__': None, '__package__': '', '__loader__': <_frozen_importlib_external.SourceFileLoader objec...
         └ <code object <module> at 0x7fa7c7265540, file "/home/hugo/.vscode/extensions/ms-python.python-2020.5.86806/pythonFiles/lib/py...

  File "/home/hugo/.vscode/extensions/ms-python.python-2020.5.86806/pythonFiles/lib/python/debugpy/no_wheels/debugpy/__main__.py", line 45, in <module>
    cli.main()
    │   └ <function main at 0x7fa7c2eed7b8>
    └ <module 'debugpy.server.cli' from '/home/hugo/.vscode/extensions/ms-python.python-2020.5.86806/pythonFiles/lib/python/debugpy...

  File "/home/hugo/.vscode/extensions/ms-python.python-2020.5.86806/pythonFiles/lib/python/debugpy/no_wheels/debugpy/../debugpy/server/cli.py", line 430, in main
    run()
    └ <function run_file at 0x7fa7c2eed598>

  File "/home/hugo/.vscode/extensions/ms-python.python-2020.5.86806/pythonFiles/lib/python/debugpy/no_wheels/debugpy/../debugpy/server/cli.py", line 267, in run_file
    runpy.run_path(options.target, run_name=compat.force_str("__main__"))
    │     │        │       │                │      └ <function force_str at 0x7fa7c5684d90>
    │     │        │       │                └ <module 'debugpy.common.compat' from '/home/hugo/.vscode/extensions/ms-python.python-2020.5.86806/pythonFiles/lib/python/debu...
    │     │        │       └ '/home/hugo/workspace/venv3/bin/goquartical'
    │     │        └ <debugpy.server.cli.Options object at 0x7fa7c2ecbf28>
    │     └ <function run_path at 0x7fa7c5c08400>
    └ <module 'runpy' from '/usr/lib/python3.6/runpy.py'>

  File "/usr/lib/python3.6/runpy.py", line 263, in run_path
    pkg_name=pkg_name, script_name=fname)
             │                     └ '/home/hugo/workspace/venv3/bin/goquartical'
             └ ''

  File "/usr/lib/python3.6/runpy.py", line 96, in _run_module_code
    mod_name, mod_spec, pkg_name, script_name)
    │         │         │         └ '/home/hugo/workspace/venv3/bin/goquartical'
    │         │         └ ''
    │         └ None
    └ '__main__'

  File "/usr/lib/python3.6/runpy.py", line 85, in _run_code
    exec(code, run_globals)
         │     └ {'__name__': '__main__', '__doc__': None, '__package__': '', '__loader__': None, '__spec__': None, '__file__': '/home/hugo/wo...
         └ <code object <module> at 0x7fa7c2ee2db0, file "/home/hugo/workspace/venv3/bin/goquartical", line 3>

> File "/home/hugo/workspace/venv3/bin/goquartical", line 11, in <module>
    load_entry_point('quartical', 'console_scripts', 'goquartical')()
    └ <function load_entry_point at 0x7fa7c3893268>

  File "/home/hugo/workspace/QuartiCal/quartical/executor.py", line 57, in execute
    model_xds_list = add_model_graph(preprocessed_xds_list, opts)
                     │               │                      └ Namespace(G_direction_dependent=False, G_freq_interval=0, G_time_interval=4, G_type='complex', _bitflag_exists=False, _bitfla...
                     │               └ [<xarray.Dataset>
                     │                 Dimensions:        (ant: 16, chan: 119, corr: 2, row: 4440, uvw: 3)
                     │                 Coordinates:
                     │                     ROWID          (row) i...
                     └ <function add_model_graph at 0x7fa77bc32598>

  File "/home/hugo/workspace/QuartiCal/quartical/data_handling/model_handler.py", line 26, in add_model_graph
    predict(ms_xds, opts) if opts._predict else [{}]*len(ms_xds)
    │       │       │        │    │                      └ [<xarray.Dataset>
    │       │       │        │    │                        Dimensions:        (ant: 16, chan: 119, corr: 2, row: 4440, uvw: 3)
    │       │       │        │    │                        Coordinates:
    │       │       │        │    │                            ROWID          (row) i...
    │       │       │        │    └ True
    │       │       │        └ Namespace(G_direction_dependent=False, G_freq_interval=0, G_time_interval=4, G_type='complex', _bitflag_exists=False, _bitfla...
    │       │       └ Namespace(G_direction_dependent=False, G_freq_interval=0, G_time_interval=4, G_type='complex', _bitflag_exists=False, _bitfla...
    │       └ [<xarray.Dataset>
    │         Dimensions:        (ant: 16, chan: 119, corr: 2, row: 4440, uvw: 3)
    │         Coordinates:
    │             ROWID          (row) i...
    └ <function predict at 0x7fa76a98ad08>

  File "/home/hugo/workspace/QuartiCal/quartical/data_handling/predict.py", line 724, in predict
    for stype in group_sources.keys()]
                 │             └ <method 'keys' of 'dict' objects>
                 └ {'gauss': Gauss(radec=dask.array<array, shape=(1364, 2), dtype=float64, chunksize=(10, 2), chunktype=numpy.ndarray>, stokes=d...

  File "/home/hugo/workspace/QuartiCal/quartical/data_handling/predict.py", line 724, in <listcomp>
    for stype in group_sources.keys()]
        │        │             └ <method 'keys' of 'dict' objects>
        │        └ {'gauss': Gauss(radec=dask.array<array, shape=(1364, 2), dtype=float64, chunksize=(10, 2), chunktype=numpy.ndarray>, stokes=d...
        └ 'gauss'

  File "/home/hugo/workspace/QuartiCal/quartical/data_handling/predict.py", line 661, in vis_factory
    dde, jones, dde, die, None, die)
    │    │      │    │          └ dask.array<broadcast_to, shape=(61, 16, 119, 2, 2), dtype=complex128, chunksize=(47, 16, 119, 2, 2), chunktype=numpy.ndarray>
    │    │      │    └ dask.array<broadcast_to, shape=(61, 16, 119, 2, 2), dtype=complex128, chunksize=(47, 16, 119, 2, 2), chunktype=numpy.ndarray>
    │    │      └ None
    │    └ dask.array<chunk_einsum, shape=(1364, 4440, 119, 2), dtype=complex128, chunksize=(10, 3480, 119, 2), chunktype=numpy.ndarray>
    └ None

  File "/home/hugo/workspace/codex-africanus/africanus/rime/dask_predict.py", line 362, in predict_vis
    die1_jones, base_vis, die2_jones)
    │           │         └ dask.array<broadcast_to, shape=(61, 16, 119, 2, 2), dtype=complex128, chunksize=(47, 16, 119, 2, 2), chunktype=numpy.ndarray>
    │           └ None
    └ dask.array<broadcast_to, shape=(61, 16, 119, 2, 2), dtype=complex128, chunksize=(47, 16, 119, 2, 2), chunktype=numpy.ndarray>

  File "/home/hugo/workspace/codex-africanus/africanus/rime/predict.py", line 419, in predict_checks
    raise ValueError("One of the following pre-conditions is broken "

ValueError: One of the following pre-conditions is broken (missing values are ignored):
dde_jones{1,2}.ndim == source_coh.ndim + 1
dde_jones{1,2}.ndim == base_vis.ndim + 2
dde_jones{1,2}.ndim == die_jones{1,2}.ndim + 1

One direction independent gain is being solved for (works when specifying MODEL_DATA)

Only have_coh and have_dies are true

Expected sizes are not the same, which is why the check is failing

expected_sizes
[[5, 4, 3, 4], [6, 5, 4, 5]]

Now I'm quite stuck as I'm not sure what this routine is actually supposed to be checking.

[bennahugo]

I was hoping to write the fft routines for lsms first to get a feel for how things stick together in the package, placing components onto a uniformly sampled grid enough to encompass the entire lsm field, and then degrid into multiple facets with my geometric routines, but will need help with this before continuing

[JSKenyon]

Can you paste your command line please?

[bennahugo]

    "configurations": [
        {
            "name": "Quartical launch (python)",
            "type": "python",
            "request": "launch",
            "program": "/home/hugo/workspace/venv3/bin/goquartical",
            "args": ["--input-ms-name", "/home/hugo/workspace/msdir/1491291289.1ghz.1.1ghz.4hrs.ms",
                     "--input-ms-column", "DATA",
                     "--input-ms-weight-column", "UNITY",
                     "--input-ms-time-chunk", "300s",
                     "--input-ms-freq-chunk", "0",
                     "--output-visibility-product", "predicted_model",
                     "--output-column", "MODEL_OUT",
                     "--solver-gain-terms", "G",
                     "--parallel-nthread", "8",
                     "--parallel-scheduler", "threads",
                     "--G-type", "complex",
                     "--G-time-interval", "4",
                     "--G-freq-interval", "0",
                     //"--input-model-recipe", "/home/hugo/workspace/DEEP2_tagged.lsm.html+~/home/hugo/workspace/DEEP2_tagged.lsm.html@dE:/home/hugo/workspace/DEEP2_tagged.lsm.html@dE",
                     "--input-model-recipe", "/home/hugo/workspace/DEEP2_tagged.lsm.html",
                     "--flags-write-out", "False"],
            "console": "integratedTerminal",
            "env": {"PYTHONPATH": "${workspaceFolder}",
                    "PATH":"${workspaceFolder}/../venv3/bin"},
            "pythonPath": "/home/hugo/workspace/venv3/bin/python",
            "cwd": "${workspaceFolder}/.."
        }

[bennahugo]

Data does not have cross correlations btw. Not sure if this is the problem though

[JSKenyon]

Was about to say, I think that might be it. How big is the dataset? If you can make it available to me I can take a look on Monday.

[bennahugo]

About 1.1G. I can put it on the ftp

[bennahugo]

Note: there aren't any delta scale components in the lsm though

[bennahugo]

Confirmed. You can replicate the bug for a dataset generated using

simms -dir "J2000,04h13m20.0s,-80d00m00s" -T meerkat -dt 16 -st 0.5 -sl 0.16667 -nc 512 -f0 856MHz -df 1044.9215kHz -pl XX YY -n mk64.0.5hr.10min.scan.16s.1045Mhz.deep2.ms

but not for a full polarimetric dataset generated using

simms -dir "J2000,04h13m20.0s,-80d00m00s" -T meerkat -dt 16 -st 0.5 -sl 0.16667 -nc 512 -f0 856MHz -df 1044.9215kHz -pl XX XY YX YY -n mk64.0.5hr.10min.scan.16s.1045Mhz.deep2.ms

[bennahugo]

I will continue with full polarimetric datasets for now as I'm only interested in building the prediction system.

[bennahugo]

I think one has to first pad the data up to full resolution to be able to apply jones matricies? At least I suspect that things are going wrong in https://github.com/JSKenyon/QuartiCal/blob/master/quartical/data_handling/predict.py#L647 As far as I can see it returns a shape of dimension 5 when 4 correlations are available and a shape 4 when 2 correlations are used

[JSKenyon]

Digging into this a little - this is currently failing due to P-jones construction, which is likely my fault. Shouldn't be too difficult to fix. Will keep you updated.

[bennahugo]

It is a bit weird though as I said - even for a diagonal matrix one cannot simply express it as two elements: what if one is calibrating the leakages or crosshand phase using only the crosshands (ie. the ms only has XY and YX - this is quite possible to have), or even more generic one only has RR / XX / LL / YY. This needs to be put in the correct place in a 2x2 matrix in order to apply the correct corruption? I think that this really does need to be in a 2x2 form for further gains to be applied correctly.

[JSKenyon]

Not necessarily. I already work with a flat correlation axis during calibration. The important point is how to consistently predict the correct values given differing numbers of correlations in the MS. The user may want to predict/calibrate with fewer correlations than are in the MS. This is sort of non-trivial in the predict layer. The easiest approach would be to always predict four correlations and then select out, but I worry that it is wasteful.

[bennahugo]

Ok, as long as those elements make it into the correct indexes in that flat array. To be honest though the way codex handles 4 correlation predict doesn't really make sense as it stands. The Q and U terms really needs a much higher order polynomial / RM column to actually predict as it is far more structured than just a spectral index. So one could restrict this to I only for now and broadcast - no functionality will be lost in the process.

[bennahugo]

I'm only planning to support stokes I for now in the fft predict and broadcast it into the correlations matrix.

[JSKenyon]

I think I will wait for @sjperkins to be available and then we can discuss it in detail. My understanding of things like feed rotation fall down a bit when we look at two correlation data.

[JSKenyon]

As a more general comment, I think I need to make sure that two correlation data is consistently handled throughout the predict. I think it is unreliable at the moment, largely as a result of me having copied a slew of code from codex, without necessarily adapting it perfectly.

[bennahugo]

To add to that one would predict it as follows:

V = K(HA) \outer (I(nu) . (|P| . e^{i . 2(RM.\lambda^{2} + p)}))

where p is the EVPA (0.5 arctan(U/Q)) and RM is given in rad / m^2 and P is the fractional polarization (-1 < |P| < 1, P := U + iQ). K is the scalar Muller term describing the baseline delay in the direction of the source. This entire term is time-invariant, so one does not actually use the DFT to predict I, Q, U and V. One only predicts I and then infer Q, U and V from that. V is also a time and frequency invariant fraction. One never does a 4 correlation predict from a fitted catalog!

The only time when it becomes time-variant is when one applies the sky rotation P term onto it.

One can additionally add a Taylor polynomial (NOT log polynomal) to describe any depolarization of the source over really wide bandwidths.

The DFT is currently the most expensive thing in the entire calibration, so you can cut your runtime by 4x @JSKenyon

[bennahugo]

Ok Isee the K term is computed once in cubical. However the note on the rest of the terms remain valid to describe the full polarimetric response

[JSKenyon]

Note that this should definitely be less broken than it once was, but I am going to leave this open until such time as we figure out precisely what we want to do in the polarised case. Thanks again @bennahugo for all your contributions here.