Discuss SPI fitter functionality

[landmanbester]

Just transferring the discussion from africanus.

Currently, the use cases are as follows:

Make power beam pattern

This is such a ubiquitous operation, this should be a standalone script. Arguments:

• image, cube or MFS image, this is used as a donor for the FITS header, and sets up the L/M/FREQ grid (and STOKES, when we get to polarization eventually)

• --ms MS --field FIELD, as per current script, used to get range of PA angles.

• --sparsify-time N, as per current script, takes every Nth timeslot

• -o OUTPUT_FILENAME, output file (should be IMAGE.pb.fits by default)

• --beam-model FILENAME, base filename of beam model

For simplicity, you could even decree that this functionality only lives in the standalone script (i.e. other scripts would just read a PB cube produced by this script, and not try to make their own.)

Convolve and/or correct single image

Take input image (or cube), and convolve it to a specified resolution and/or apply a primary beam correction (either, or both).

• image, image to work with

• --psf MAJ MIN PA, desired PSF parameters (prefer to use PSF over "beam" whatever Steve Gull says, because otherwise there's too many "beams" in this pot).

• --circ-psf MAJ, alternative to --psf, use circular PSF of the given size (called --circularize_beam now, but I think this is better).

• -o OUTPUT_FILENAME, output file

Read image, read its BMAJ/BMIN/BPA, work out what additional convolution to apply to make the resolution match the specified MAJ/MIN/PA, apply the convolution. Write the new PSF size to the FITS header. If the image has no BMAJ/BMIN/BPA, assume it is a model image and just apply a convolution by specified PSF directly.

• --pb PBIMAGE, if specified, apply a power beam correction from the given beam file (which must match the image size).

• --pb-min 0.05 set output to 0 where power beam drops below this value.

Spectral index fitter

• model [residual] model & residual cubes

• --psf or --circ-psf optional. If specified, convolves to specified resolution (otherwise use resolution of lowest channel?)

• --pb PBIMAGE.fits optional. If specified, applies power beam correction using the specified beam cube (and apply --pb-min as specified), and also applies it to the I0 map.

Other options as per current script. Except I don't like the use of --output to specify the set of output images (it's more conventionally used to specify an output filename). How about

• --products which images to write (currently known as --output)

• -o/--output base output filename (default should be to take model and strip the .fits extension)

[o-smirnov]

Suggest making -st 10 the default for the power beam maker. Every timeslot takes too long, and I can't imagine that sort of accuracy is necessary. (As it is, took me 10m with -st 10). Also, the option needs to be documented a little better.

[o-smirnov]

simple_spi_fitter wants pypocketfft. I assume I can just install from the repo https://gitlab.mpcdf.mpg.de/mtr/pypocketfft? (Can we persuade Martin to release it on PyPI and just have it installed via setup.py?)

[landmanbester]

Wait, those exact functions are in ducc0 now, I'll change the imports quickly

[landmanbester]

All done in the spi_fitter branch. Can you please give it a go?

[o-smirnov]

Now getting:

    from ducc0 import r2c, c2r
ImportError: cannot import name 'r2c'

...despite

$ pip install ducc0
Requirement already satisfied: ducc0 in /home/oms/.venv/cc/lib/python3.6/site-packages (0.7.0)

[landmanbester]

Too mush haste too little speed. Try again

[o-smirnov]

Yep, works now.

[o-smirnov]

@landmanbester, could you take a look at /net/young//home/oms/projects/Shapley/selfcal-shap1/a3562 when you have a chance?

Running it as

simple_spi_fitter.py -model a3562-model.fits -residual a3562-residual.fits -o a3562-spi -bm a3562-beam.fits -th 0.5

Even with the threshold that low, it's being very picky about what it fits:

restored image	I0 image

We obviously want to get on top of that fluffy stuff...

[o-smirnov]

Aha, I think I see the problem. The convolved model is way too faint (compared to the corresponding restored image), and the clean_psf image has a peak of 3e-3. This convolution function needs to have a peak of 1, right? I guess the fix for #36 broke this one...

[landmanbester]

This convolution function needs to have a peak of 1, right?

Yep, I just pushed a fix.

We obviously want to get on top of that fluffy stuff...

If you have a lot of flux left in the residual that could be problematic. I wonder if it is worth/wise to also convolve your residuals to a common resolution and add them back in. That means you will have some contribution from psf sidelobes but at least you are not missing the flux altogether. It really is time to do a realistic simulation with extended sources to see what works best. I have those resolve images of Cygnus A from 2-13 GHz that could be used for this. Let me see if I can set something up. In the meantime, is it worth adding an option to convolve and add in the residuals?

[o-smirnov]

In the meantime, is it worth adding an option to convolve and add in the residuals?

If it's relatively easy, yes please. Sounds like the thing to do for a lot of faint fluffy stuff.

[o-smirnov]

Ah, an embarrassment of riches now, component-wise, but...

Fitting 273393 components
Warning - max iterations exceeded for component  272
Traceback (most recent call last):
  File "/home/oms/.venv/cc/bin/simple_spi_fitter.py", line 7, in <module>
    exec(compile(f.read(), __file__, 'exec'))
  File "/home/oms/projects/pfb-clean/scripts/simple_spi_fitter.py", line 368, in <module>
    main(args)
  File "/home/oms/projects/pfb-clean/scripts/simple_spi_fitter.py", line 283, in main
    np.float64(ref_freq)).compute()
  File "/home/oms/.venv/cc/lib/python3.6/site-packages/dask/base.py", line 279, in compute
    (result,) = compute(self, traverse=False, **kwargs)
  File "/home/oms/.venv/cc/lib/python3.6/site-packages/dask/base.py", line 567, in compute
    results = schedule(dsk, keys, **kwargs)
  File "/home/oms/.venv/cc/lib/python3.6/site-packages/dask/threaded.py", line 84, in get
    **kwargs
  File "/home/oms/.venv/cc/lib/python3.6/site-packages/dask/local.py", line 486, in get_async
    raise_exception(exc, tb)
  File "/home/oms/.venv/cc/lib/python3.6/site-packages/dask/local.py", line 316, in reraise
    raise exc
  File "/home/oms/.venv/cc/lib/python3.6/site-packages/dask/local.py", line 222, in execute_task
    result = _execute_task(task, data)
  File "/home/oms/.venv/cc/lib/python3.6/site-packages/dask/core.py", line 121, in _execute_task
    return func(*(_execute_task(a, cache) for a in args))
  File "/home/oms/.venv/cc/lib/python3.6/site-packages/dask/optimization.py", line 963, in __call__
    return core.get(self.dsk, self.outkey, dict(zip(self.inkeys, args)))
  File "/home/oms/.venv/cc/lib/python3.6/site-packages/dask/core.py", line 151, in get
    result = _execute_task(task, cache)
  File "/home/oms/.venv/cc/lib/python3.6/site-packages/dask/core.py", line 121, in _execute_task
    return func(*(_execute_task(a, cache) for a in args))
  File "/home/oms/.venv/cc/lib/python3.6/site-packages/africanus/model/spi/dask.py", line 27, in _fit_spi_components_wrapper
    maxiter=maxiter)
  File "/home/oms/.venv/cc/lib/python3.6/site-packages/africanus/model/spi/component_spi.py", line 71, in fit_spi_components
    jac, ncomps, nfreqs, tol, maxiter)
ZeroDivisionError: division by zero
Warning - max iterations exceeded for component  802

[landmanbester]

Ah man, an annoyance I thought I got rid of. It happens when the GN optimisation goes awry. I've been meaning to switch to LM instead but maybe it's worth checking why it goes wrong in the first place. What settings did you run it with? I hope you changed -th to something more sensible?

[o-smirnov]

I had it at -th 5. At -th 10 I get 145346 components and no errors. So yeah, it's probably just a low-SNR component going awry. Might be better to just skip it then and issue a warning at the end, something like "N components failed to fit, you might want to set your thresholds higher".

[landmanbester]

Yeah, and maybe set the default to 10. 5 is a bit low for reliable spi's (if there even is such a thing)

[o-smirnov]

I'm not getting any of the fluffy stuff at th=10 though. So I think a check for the fitting error, and convolving the residuals in, might be the only way to get something. Even if it's super unreliable.

[landmanbester]

Ok, I'm working on it. An easy hack to add the residuals back in for now is to convolve the restored image to a common resolution and pass the --dont-convolve flag to the spi fitter

[landmanbester]

I'm not getting any of the fluffy stuff at th=10 though

Do you see the fluffy stuff in the convolved model though? Only written out of m is one of the products

[o-smirnov]

Some of it is in the model, but it's probably falling below the threshold:	conv model	spi

[landmanbester]

So I think a check for the fitting error, and convolving the residuals in, might be the only way to get something. Even if it's super unreliable.

Checking for the fitting error is more difficult than it sounds because neither try/except nor if statements are supported inside compiled numba code where the check has to happen. I've tried to hack around it by simply hard coding a minimum value for the determinant here but I'm not sure if that is going to result in nonsense or not. I'll test it presently but if you want to give it a go so long you can checkout my fix_div0_spi branch.

A better fix for this would be to regularise the inversion towards some mean (-0.7?) and then use the approximate posterior uncertainties as a reliability indicator. I am busy working on that too but for some reason my uncertainty estimates are way off at the moment and I still can't figure out why.

In any case, I have to get this working for @svw26's data so expect some changes to the interface.

[o-smirnov]

OK well the good news is that I can now fit an insane number of components (-th 2), and it even starts picking out the fluffy stuff, and the results aren't even crazy:

Looking forward to trying this with convolved residuals in, might even out the "bubbly" texture in the extended stuff a bit.

The spi_err map has a scaling problem though -- it's all in the 1e+5 -- 1e+8 range. Qualitatively ok in that it seems to follow the SNR, just the scaling is weird.

[landmanbester]

Cool, thanks for testing @o-smirnov. I'll have the convolved residuals in there shortly.

The spi_err map has a scaling problem though -- it's all in the 1e+5 -- 1e+8 range. Qualitatively ok in that it seems to follow the SNR, just the scaling is weird.

Yeah, busy trying to track this down. What does the I0err map look like? If we use the actual wsums (which are written out to the fits header by wsclean but lost during stacking the cube) this can probably act as a poor man's uncertainty estimate...

[o-smirnov]

The I0err is not as crazy (0 to 5 range), but still seems high. An error of 2 Jy on the brightest (6 mJy, i.e. SNR 1000) component can't be right.

[landmanbester]

The problem with I0var and alphavar stemmed from not providing absolute variances on the observed per band intensity estimates. The workaround is to insist that the reduced chi square equates to one which adds a simple scaling to the estimates. I've added the scaling to the scrips in africanus so the lower bounds on the variances in alpha and I0 should be approximately "correct" now (use latest africanus master branch). I've also refactored things a bit, convolved residuals coming soon

[landmanbester]

Alright, if a residual cube is provided to simple_spi_fitter.py it will convolve it to the same resolution as the convolved model before fitting spectral indices. I have refactored the code a little and tested all three the power_beam_maker, image_convolver and simple_spi_fitter locally. They all seem to run through but I haven't tested how reasonable the results are. I would be grateful if you could give this a go with your data @o-smirnov.

Some TODOs:

• work directly on wsclean output (i.e. bypass stacking the cube),
• work directly with eidos coefficients instead of interpolating a fits file (hoping this functionality will make it into africanus soon),
• set up realistic reference simulation (maybe using Cygnus A?) and implement and end to end test.

Any other (reasonable) requests?

[o-smirnov]

OK I updated and rerun, but I got numerically identical results for the spi image, which seems unlikely. Are you sure the residuals are being included in the fit by default?

[landmanbester]

Yep, the default is to add the convolved residual in if the residual is passed in and it finds the beampars in the header. Did you see "Convolved residuals added to convolved model" or "Can't find Gausspars in residual header, unable to add residuals back in" in the log?

[o-smirnov]

Not that I can see: https://pastebin.com/y7M35EqW

My codex is up-to-date with lbester/fix_div0_spi, and pfb-clean with spi_fitter.

[landmanbester]

Are you sure you sure you updated the spi_fitter branch? The log seems out of date. In particular it won't print "Doing FFT's" anymore. The changes to africanus have been merged into the latest master

[o-smirnov]

Ah, I know what it is. I need to pip install -e again. The "editable" install symlinks everything except binary scripts -- those get copied over.

[o-smirnov]

Hmmm, that didn't go well... the SPI image seems upside down now, I'm getting way more components for the same settings, and the values are all over the shop:

before	after

[landmanbester]

Sigh, I need test cases! Let me see what's going wrong

[o-smirnov]

Might be as simple as your model coordinates being inverted? The structure looks flipped in both X and Y.

[landmanbester]

I'm not sure why but for some reason the spatial fits axes come out flipped when the data is 4D (i.e. has a Stokes axis). This is probably a bug in the way I handle fits files but I don't really want to change it right now because all other scripts use this convention. Will have to figure out what's going wrong at a later stage.

As for the crazy spectra, I guess this would happen if the threshold or pb cutoff is set too low. What values were you running with?

[o-smirnov]

I have a -th 2 which I guess is pretty low, but I had the same in the "before" case and the result wasn't as wild.

Maybe adding in the residuals makes things worse for weak components. is there a command-line option to disable that again? (Because if I don't supply residuals, how do you estimate thresholds?)

[landmanbester]

There is now. It will only add them in if you pass the -acr flag. If no residuals are passed in then I use the maxDR option to set thresholds but getting this from the residual is the better way to go by far.

I am just realising that dividing by the clean beam when residuals have been added in is probably not ideal. It will have a similar effect on the noise as "correcting" your data. We should actually put the beam in the forward model and use that to fit spi's. This way the uncertainty estimates should reflect the SNR as well. Shouldn't be too difficult to do, will try it and see if it improves things

[o-smirnov]

OK thanks. Sadly, I can confirm now that adding the residuals has made it worse. Without -acr I get the same image as before. With -acr, things get funky. Not sure whether there's remaining problems, or because this was just a bad idea to begin with, but:

no -acr	-acr

...check out the mottled dark features across the compact positive-spi source at 7 o'clock. Zooming in:

spi	convolved_model	convolved_residual

I'm very suspicious of those high-freq ripples across the residual. It shouldn't be there in a convolved image, right? Something with the padding?

[o-smirnov]

I also wonder -- stop me if this is nonsense -- should we try a mode where we detect components in the model only, and fit in model+residuals?

[landmanbester]

I'm very suspicious of those high-freq ripples across the residual. It shouldn't be there in a convolved image, right? Something with the padding?

Very suspicious. The convolved residuals should also be written out so you can add them together manually and see what the spectrum looks like. It is very possible that something is going wrong. Are those components anywhere near the null of the beam?

I also wonder -- stop me if this is nonsense -- should we try a mode where we detect components in the model only, and fit in model+residuals?

I don't know if you recall but that is what we did initially. The problem comes in when you want to make a map because the individual spectra are not at all physical. We ended up weighting the derived components by their fluxes before convolving them. That gave pretty good results but we were at a loss to figure out what to do with the uncertainty maps. I think the way to go is to try with the beam in the forward model (as apposed to dividing by it). Shouldn't be a major change. Do you need ithese maps urgently?

Edit - Ah, you are looking at the convolved residual. Let me check something quickly

[landmanbester]

Yeah, could be aliasing. The residual does not go to zero towards the edges (although the convolution kernel is small so I thought that wouldn't be a problem) . Maybe try with -pf 2. In the meantime, let me write some tests to make sure things are working as they should

[o-smirnov]

Are those components anywhere near the null of the beam?

On the contrary -- quite central.

[landmanbester]

After a discussion with @IanHeywood I realised that the SPI fitter didn't convolve the residuals correctly. I was normalising the convolution kernels to have a volume of one so that flux is conserved when you convolve them to a common resolution. This will artificially lower the the flux density at higher frequencies resulting in biased spectra. Normalising the kernels to have a unity peak automatically takes care of the flux scaling when adjusting the resolution. I have pushed a fix to the refactor_minors branch. I have also added a test which ensures that this is the correct thing to do https://github.com/ratt-ru/pfb-clean/blob/refactor_minors/pfb/test/test_convolve2gaussres.py. I'm not sure how much this affects your spi maps @o-smirnov but might be worth checking

[o-smirnov]

Hmmm, quite possible -- Tiziana was reporting some differences w.r.t. fitting in AIPS -- let me try again!

[landmanbester]

Now discussed! I have transferred this functionality to spimple which is also available on pypi now. Let's continue the discussion there if needs be