very slow convergence

[o-smirnov]

I keep breaking it... I've been noticing things getting slower and slower as I get deeper into this data, but this latest one is beyond ridiculous. As you can see from the log below, it's taking on the order of several hours -- on a GPU -- to do a single 2Kx2K image, with a 1Kx1K subregion (and a mask, which should restrict things even further). As usual, the MS is in /net/aretha/home/oms/VLA-CygA/MS-Jon (NB: it's the "HI" MS this time), and the mask is one directory up from that.

$ grep Parameters reduction-21Jul15-HI-224chan-chan16-lsmthr0.25-spi-moresane3-rob-2/log-CYG-S-HI.txt 
[07/21/2015] [04:44:22] [INFO]: Parameters:
[07/21/2015] [08:16:15] [INFO]: Parameters:
[07/21/2015] [11:17:10] [INFO]: Parameters:
[07/22/2015] [01:32:14] [INFO]: Parameters:
[07/22/2015] [01:54:51] [INFO]: Parameters:
[07/22/2015] [02:21:02] [INFO]: Parameters:
[07/22/2015] [02:50:27] [INFO]: Parameters:
[07/22/2015] [06:18:49] [INFO]: Parameters:
[07/22/2015] [11:38:23] [INFO]: Parameters:

Here's the command:

/home/oms/src/wsclean-code/wsclean/build/wsclean -niter 3 -mgain 0.85 -scale 5.555
55555556e-05 -fitsmask mask-C-HI-2048-exp.fits -superweight 4 -weight briggs -2.00 -moresane-arg '-sbr 1024 --negcomp --
mask mask-C-HI-2048-exp.fits -sl 5' -channelrange 16 240 -datacolumn CORRECTED_DATA -name reduction-21Jul15-HI-224chan-c
han16-lsmthr0.25-spi-moresane3-rob-2/plots-CYG-S-HI-128-224/CYG-S-HI-s2.moresane -field 2 -pol IQUV -gain 0.05 -moresane
-ext /home/oms/PyMORESANE/pymoresane/bin/runsane -threshold 0 -nosmallinversion -size 2048 2048 -channelsout 14 CYG-S-HI.MS

I suspect it may be a convergence problem induced by us enabling negative components -- and also this being very "noisy" data due to relatively bad calibration. I'll try it without negcomp, and maybe with a higher sigma level...

[o-smirnov]

A further oddity is that, despite the -sl 5 setting above, it seems to carry on modelling well into the "noise". Look at these residuals:

With noisy data like that, I would like to stop deconvolution well before it starts trying to pick up on the artefacts (analogy being using a high stopping threshold with normal CLEAN). I thought the --sigmalevel option would do that, but apparently not. Any suggestions?

[JSKenyon]

This is a problem with the masking procedure. At the moment, the mask zeros out all the values which do not fall inside it before estimating the noise. This is done so that the wavelet scales are only constructed with the data in the mask. The problem is that this means that the noise is incorrectly estimated and as a result you get this overcleaning behaviour. I will think about the problem and try and devise a solution. I have an idea, but it may be computationally expensive.

[JSKenyon]

To clarify, it is cleaning to -sl 5, but the sl in the masked region is obviously much less than the sl outside the masked region.

[JSKenyon]

A cheat which may work in the interim is to set -sl high. An appropriate value could be approximated by estimating the noise outside the mask and comparing it to the noise inside the mask and picking a -sl value to make them comparable.

[o-smirnov]

Yeah I suspected something like that... well in this case the correct thing to do is simply to estimate noise by going outside the mask. In fact I think this should be the default behaviour. After all, the whole point of the mask is to tell the algorithm where the emission is (and isn't). Thus the non-masked regions are exactly where we should be estimating the noise.

[JSKenyon]

The issue comes in estimating the noise in the wavelet scales themselves, as then the mask is a different size in each scale. My idea for dealing with this is to decompose the mask and create a three dimensional array of masks which you apply to the wavelet decomposition. This is the potential computation expense I was referring to.

[JSKenyon]

I am having difficulty getting the dirty image and PSF. I cannot run it on Aretha and Jake no longer has space for me to copy the MS into. Could you possibly put the dirty image and PSF somewhere? Think it will be easier than trying to find space for the MS.

[JSKenyon]

OK, this masking thing is less simple that I would have thought. Decomposing the mask doesn't really work as you end up with all sorts of nonsense at the large scales. The same is true for simply dilating the mask by the filter width at each scale. The current system works better, but needs tweaking. I think the best approach will be to estimate the noise from a full decomposition but work only on the masked data. This shouldn't be too difficult, but may require a bit of fiddling.

[o-smirnov]

Yep, noise from full image sounds like the best plan. I shall desist from using masks for now. Let me know if/when you've got a quick fix.

[JSKenyon]

I have written something which will hopefully do what you need. Could you possibly give me the dirty image, PSF and mask in question - I don't have particularly good test data (no space on jake to copy the MS into).

[JSKenyon]

Should be fixed in the latest release - I am fairly confident it is working. Please let me know if there are any problems with it.