Creating the Cyg X-3, focus=1903MHz LSS reference image.

The basic problem: at 1903, the PB is extra wide. The bright sources
around nominal half-power are even brighter. And it seems that
lwimager/maximum-entropy doesn't deconvolve them right, giving really
nasty corrugations in the maxen model.

Demonstration: rephase the dataset to be centered on the bright
sources.  Image with lwimager/maxen. Residuals near the bright sources
are much better, as is the model.

mfsimager wprojplanes=128 npix=2048 cellsize=10arcsec operation=entropy \
  sigma=1mJy prior=../miscdata/memprior-pb2mJy.img srls.ms image=t.rph1 \
  phasecenter='j2000, 20h32m31.2s, 40d16m20s'
mfsimager wprojplanes=128 npix=2048 cellsize=10arcsec operation=entropy \
  sigma=1mJy prior=../miscdata/memprior-pb2mJy.img srls.ms image=t.rph2 \
  phasecenter='j2000, 20h28m06.6s, 40d51m40s'
mfsimager wprojplanes=128 npix=2048 cellsize=10arcsec operation=entropy \
  sigma=1mJy prior=../miscdata/memprior-pb2mJy.img srls.ms image=t.rph3 \
  phasecenter='j2000, 20h33m55.5s, 41d22m46s'

Tedious, dissatisfying solution: use these rephased images and
deconvolution solutions to subtract out the bright guys, use maxen to get
the rest of the flux, and merge everything back together.

Problems with the solution: most of it is just the hassle, and being sad
that things don't just work right out of the box, but I can't find a 
good way to recombine the various model components. I ended up doing some
manual image processing as will be described shortly.

Step 1: Extract the relevant sources out of the rephased images.

msimsub 958 949 1124 1095 t.rph1.model t.rph1.model.foc
msimsub 993 966 1106 1058 t.rph2.model t.rph2.model.foc
msimsub 993 966 1106 1058 t.rph3.model t.rph3.model.foc

Step 2: Experiments indicated that these little postage stamps were
being treated as being surrounded by zeros, so if you just subtract
the above models from the data you get a big negative box for each
stamp. Therefore you need to subtract off the "typical" background
flux value from each stamp. Both mean and median aren't good
estimators since the whole point is that the stamp is filled with a
bright extended source. I plotted histograms and chose values by eye.
I then modified the images by loading them into ipython and futzing.
Something like

$ cp -r t.rph1.model.foc{,.mnsub}
$ ipython
> import numpy as np, astimage, omega as om, ndshow
> i = astimage.open ('t.rph1.model.foc.mnsub', 'rw').simple ()
> d = i.read ()
> om.quickHist (d.flatten (), range=(0, 0.006), bins=40).show ()
> # choose a value
> d -= typical_bg_value
> i.write (d)
> i.close ()
> # repeat 2 more times

Step 3: Subtract the postage stamps from the UV data. This requires
some craftiness but fortunately mswpftsub lets us do the right thing.
Points are:

- by default, mswpftsub resets the model data in the ft step
- the mswpftsub step is cumulative on existing corrected_data

The initial dataset has the model filled with the best point source
guess (BTW, this was all based on setupreclss data, since the
restorelss results were crappy because the whole point is that the LSS
was looking ugly) and the corrected_data being the raw data with the
point source models subtracted. Therefore:

mswpftsub srls.ms t.rph1.model.foc.mnsub && 
mswpftsub -i srls.ms t.rph2.model.foc.mnsub &&
mswpftsub -i -s srls.ms t.rph3.model.foc.mnsub

First we reset the model data to the first stamp, we incrementally
(-i) modify it twice, then we subtract from the preexisting
point-source-subtracted data.

Step 4: Image the rest.

mfsimager wprojplanes=128 npix=2048 cellsize=10arcsec operation=entropy \
  sigma=1mJy prior=../miscdata/memprior-pb2mJy.img srls.ms image=t.srall

Step 5: Regrid stamps to match up with the Step 4 model. (They were generated
with the alternate phase centers so the coordinate systems aren't the same.)

$ casapy
> imagename = 't.rph1.model.foc.mnsub'
> template = 't.srall.model'
> output = imagename + '.rg'
> axes = [0, 1] # things didn't work without these last two
> shape = [2048, 2048, 1, 1] # for some reason
> # repeat two more times

Step 6: Merge the models.

I was hoping that I could add things together and rerun maxen with a
good prior, or something, but it seemed to just revert to the results
of a straight maxen without any of this monkey business. Too bad but
not surprising I guess -- I'm not sure how you'd construct a
wfhogbom-type maxen algorithm that periodically went back to the
visibilities and reimaged. (I spent some time looking at the
casarest/synthesis source before deciding there was no way I could
make that happen in any noninsane amount of time.)

So, manual photoshopping in ipython.

$ cp -r t.srall.model q.model
ipython:
import astimage, numpy as np, ndshow, omega as om
ifull = astimage.open ('q.model', 'rw').simple ()
d = ifull.read ()
nd = d.copy ()
# load stamp
i1 = astimage.open ('t.rph1.model.foc.mnsub.rg', 'r').simple ()
d1 = i1.read ()
# figure out where it's defined from its mask
w = np.where (~d1.mask)
w[0].min (), w[0].max (), w[1].min (), w[1].max ()
s1 = (slice (w[0].min (), w[0].max ()+1), slice (w[1].min (), w[1].max ()+1))
# some have slight rotations, unmask, will be all zeros
d1.mask[s1] = 0
# weights to smoothly fade in the stamp model to the existing model
# ramping up 0.1 ... 0.9 along all edges
w1 = np.ones ((s1[0].stop - s1[0].start, s1[1].stop - s1[1].start))
for i in xrange (9): w1[i] = np.minimum (w1[i], 0.1 * (i + 1))
for i in xrange (9): w1[-i-1] = np.minimum (w1[-i-1], 0.1 * (i + 1))
for i in xrange (9): w1[:,i] = np.minimum (w1[:,i], 0.1 * (i + 1))
for i in xrange (9): w1[:,-i-1] = np.minimum (w1[:,-i-1], 0.1 * (i + 1))
# add the stamp model
nd[s1] = d[s1] + w1 * (d1[s1] + 0.00025)
# ... repeat two more times ...
nd[s2] = d[s2] + w2 * (d2[s2] + 0)
nd[s3] = d[s3] + w3 * (d3[s3] + 0)
# all done
ifull.write (nd)
ifull.close ()

Blinking with the plain maxen model shows lots of improvement.  Test
epimaging (after replacing fieldinfo/cygx3/lss-1903.model) shows that
things look OK (you can definitely see where the splicing has happened
in the model). Test restorelss looks WAY better, so that's good.

Unfortunately executing another iteration of epoch imaging won't make
this process any less necessary. It looks like the maxen algorithm can
be seeded with a model and maybe that's what makes the difference, so
we could at least iterate. (I tried feeding maxen the naively-combined
stuff with prior=. The semantics of model= to lwimager are very unclear
to me.)