Overview
CVS path: LOFAR/Timba/WH/contrib/OMS/fringes
Scripts for simulating data and playing with fringe fitting.
Contents
Main TDL scripts:
simul.py: simulator TDL script. This uses a source model from models.py to simulate a "clean" sky (i.e. no instrumental effects) and write it to an MS column. Use this to fill the MODEL_DATA column of an MS.
corrupt.py: this TDL script applies instrumental effects to an MS. Normal usage is to read the MODEL_DATA column, corrupt with a G Jones, and write to DATA. The included G Jones model creates a slowly varying (sine of time) phase, with a frequency slope that also evolves in time. It can also add antenna-based noise terms. The variation parameters are picked randomly per-station every time you run the script.
solve.py: this TDL script attempts to do fringe fitting on a corrupted MS. It fits the phases by a polynomial of the specified degree, then corrects the data with the inverse of the fitted G Jones. Nortmal usage is to read from the DATA column and write to the CORRECTED_DATA column.
Note that all the above scripts have a "make image" TDL job, which will run a Glish script to image the data, convert to FITS, and launch kvis to view it. It can make images in channel or MFS mode, specify this under the TDL exec menu.
Supporting scripts:
models.py: this is a collection of various source models used by simul.py to simulate a sky.
make_image.g: Glish script to run the imager. General syntax is
glish -l make_image.g <column> ms=<msname> mode=<mode>
<column> is one of DATA, MODEL_DATA, CORRECTED_DATA, or RESIDUAL. The latter option will image the difference between MODEL_DATA and CORRECTED_DATA.
<mode> is "mfs" or "channel".
<msname> is the measurement set name.
The filename for aips++ and FITS images is derived from the MS name, column and model. E.g. if the MS is called TEST.MS, column is CORRECTED_DATA and mode is "mfs", the output image will be called TEST.corrected-mfs.fits.
read_msvis_header.py: internal script called by the kernel to read an MS header and apply it to a tree.
Measurement sets:
TEST.MS: this is the CLAR MS I borrowed from Tony. Contains a VLAx10 configuration with 27 stations (so 351 baselines), 32 channels from .8 to 1.4Ghz, 8 hours at 1 minutes integration, so 480 timeslots.
Results
Fitting fringes at low SNRs
I used TEST.MS for my experiments.
1. Run simul.py with a "two point sources plus random" model. This simulates two 1 Jy point sources at fixed positions, and 25 random fainter sources (~.05-.1 Jy) all over the map, mostly point sources, with a few elliptical gaussians.
2. Run corrupt.py to corrupt the model. Specify 180 degrees for the phase excursion in both freq and time, and add .33 Jy of noise per antenna. So the effective SNR for our two bright sources is only ~3 (or perhaps slightly more? Since the combined noise per-baseline is sqrt(2) lower?)
3. Observe dirty image, we can see nothing except that there are obviously two sources, and phases are all over the place.
Here's the same image clipped at 99%:
4. Now run solve.py with the following options:
fringe fit degree: 1 (i.e. linear)
flux constraint: .5, weight 1000
source model: two point sources
solver LM factor: .001, convergence threshold .0001, max iterations: 100.
reuse solutions from previous time interval
do not reuse solutions from MEP table (i.e. solve from scratch)
other options don't really matter...
...and execute the "solve for phases and fluxes" job (load the "flux and position solutions" bookmark).
This assumes that we only know of the two point sources, and attempts a simultaneous LSQ solution (at 10-minutes intervals) for source flux and phases. Flux is assumed constant; phases are 1-degree polcs in freq and time.
The solution for the first interval is the hardest, since all the phases are completely unaligned. It can take 50+ iterations. Without the flux constraint, it can fail completely (that is, settle at some local minimum with fluxes close to 0). With the constraint in there, it eventually makes its way to the right solutions. Subsequent intervals then converge much quicker (5-8 iterations) since we already have a very close guess for the phases from the previous interval. Total runtime on birch is 4-6 minutes.
Now make the corrected image and enjoy! This is clipped at 99% to bring out the faint sources.
It is also interesting to make a residual image, and observe remaining errors. This image is clipped at 99%:
The artifacts around the bright sources range from -28mJy to +7.5mJy. Mean residual is about -20uJy. Of course some systematic errors remain, since there are 25 sources in the field that are not part of the fitted model. Despite that, this demonstrates that we can fit fringes even in extremely low-SNR situations.
Other observations
We can also fit phases with 2-degree polcs. This works fine (with barely any performance impact), but produces no appreciable difference in the maps. I'm guessing that this close to the noise level, 2nd-order phase errors don't really matter.
It would be interesting to do a different simulation: 1Jy bright sources, <=1mJy faint sources, <=1mJy noise. I expect the phase/flux fits to be much better (since the contamination from other sources is much lower), but residual phase errors should have a far greater impact on the detectability of faint sources. Here perhaps we may show a difference between 1st and 2nd-degree fringe fitting.
We can also fit phases with 2-degree polcs. This works fine (with barely any performance impact), but produces no appreciable difference in the maps. I'm guessing that this close to the noise level, 2nd-order phase errors don't really matter.
It would be interesting to do a different simulation: 1Jy bright sources, <=1mJy faint sources, <=1mJy noise. I expect the phase/flux fits to be much better (since the contamination from other sources is much lower), but residual phase errors should have a far greater impact on the detectability of faint sources. Here perhaps we may show a difference between 1st and 2nd-degree fringe fitting.
There are other models in models.py that I've played with. Some observations on those:
We can also fringe-fit on a central point source (ha ha big surprise).
Putting in a solvable spectral index really screws things up. Same for solvable positions. My guess is, chi-sq derivatives w.r.t. spectral index and position are very high compared to flux and phases, so during the initial steps the solver shoots off into the wilderness. We can probably address this once Maaijke implements simple Parm constraints -- i.e. we can then set natural bounds on position (since we always know it to within a few "pixels") and spi.
A single unaccounted-for point source at .1Jy level will screw up the solutions for the two 1Jy sources (fluxes converge to .8 rather than 1). The "two bright one faint point source" model can be used to demonstrate this. On the other hand, putting in a single extended source at the same level does not affect the solution ("two point sources plus faint extended" model). And as demonstrated above, 25 faint sources all over the map do not affect the solutions.
The reason is probably that a single point source produces a well-defined fringe, which interferes with the fringes of the two brighter sources we are modelling, and so skews the solution. The fringe from an extended source is flatter, while fringes from 25 sources will on average cancel each other out. So the moral is, as long as there's fainter flux spread all over the field, you're ok (and this is why peeling works!) On the other hand, if there's significant unaccounted-for flux localized at one (two? three?) points in the field, you have to include it in your model (and so you can't peel). (On a self-back-slapping note, this is why we need these kinds of simulations, so that we can work out calibration strategies.)
Some observations on solving in general:
Wim's recent (April 2006) changes to LSQFitter have had a definite negative impact. I had to roll back to an older aips++ build for all my experiments. In particular, with the new fitter I observed:
Much slower convergence (on the 60+ instead of 10 iteration level!!!)
Failure to stop. The new fitter is supposed to know when to stop better, but on these kinds of fits this is not the case at all. I could get it to iterate in one place with no change to chi-sq for 20-40 iterations (with epsilon=1e-4). The old solver stops much more reliably.
On the other hand, when doing ClarSimulations and fitting for beam width, flux and spectral index, the new fitter stops better than the old fitter (although it still converges slower). Go figure...
My overall conclusion is that Wim's changes to the descent algorithm are very lateral, in the sense that they may produce improvements for certain classes of chi-square surfaces, but work much worse on others. In fact, for fitting in the uv domain, there has been nothing but regression... I guess the answer is to make the descent algorithm tweakable so that we can experiment with different settings and arrive at some answers. In a perfect world, of course, we'd want Wim to play with MeqTrees directly to see how his changes impact convergence...
We need to be able to set constraints on parameters, badly. I have opened a bug on this and assigned it to Maaijke.
When solving for phases and fluxes simultaneously, we could benefit from a somewhat different fitting mode: update only phases or only fluxes at alternating iterations. This is known (I've read it somewhere but don't have a reference handy) to produce faster convergence for certain classes of problems; the phase+flux fit appears to be one of them. Can we implement this cheaply? I suppose we could tell a MeqParm to only return perturbed values at even or odd iterations, so the solver will receive two different subsets of equations at alternating iterations. This will probably interact with rank, but how? We definitely need to discuss this.