Azimuthal Averages

We noted in discussion yesterday that the aperture sum in the (nearly perfectly recovered) is lower than expected even for a gaussian fit of an airy. A Gaussian approximation to an Airy should recover 93.5% of the flux; 6.5% of the power is in the sidelobes of an Airy disk. We recover about 86-88% of the flux in the input point-source map, or about 5-7% short, despite matching the peak to within 2-3%. This discrepancy is partly - but not completely - explained by a slight negative bowl effect: the Gaussian fit has a higher total flux than the aperture sum. Using the gaussian fit value, the recovery is about 90%. A mere 3% missing flux is not too much of a problem - it could still be that the gaussian is uniformly reduced by bowl effects, or perhaps that some flux is not reincorporated into the map (at a very low level) because it is correlated enough to be PCA-cleaned.

In the figure, the top panel is the radial profile of the ds1 data (black), the input data (red), and their gaussian fits (see legend). The circles on the labeled images show the FWHM of the fitted gaussian, indicating a good fit. The apertures used for summation are NOT displayed. The sum aperture has a radius of 90" (6.25 pixels) and the "Noise Aperture" has the same area as the sum aperture.

If we redo the analysis with a smaller aperture, the recovery is better: 95% in the central 45". However, now the airy ring cannot account for any lost flux, so it still looks like there's a 5% loss intrinsic in the pipeline.

ds1-ds5 comparisons

I'm comparing simulated ds1-ds5 to real ds1-ds5 comparison tests.

In the simulated tests, I compare the recovered map after 20 iterations with 13 pca components subtracted to the input map. There are figures showing this comparison for ds1 and ds5 images in addition to one showing the comparison between ds1 and ds5. The agreement is pretty much as good as you could ask for.

These simulations are the most realistic run yet. They include a simulated atmosphere that is perfectly correlated between all bolometers excepting gaussian noise, but the relative sensitivity of the bolometers is varied.

This is what a 'real' ds1-ds5 comparison looks like. The image shown is a "cross-linked" observation of Uranus with downsampling off and on. Note that downsampling clearly and blatantly smears the source flux.

The same image with "beam location correction" looks no better.

The problem is essentially the same with the individual scan directions:

What is causing this difference?

• higher-order corrections to the atmosphere calculation?
  
• inadequate sampling of the model?
  
• "pointing" offsets between the model and the data (note that these are NOT pointing offsets, but they may be "distortion map" offsets)?
  
• Other?
  

Examining the weights and scales for two individual (real) observations, ds1 followed by ds5, is not particularly telling; there is one additional outlier bolometer flagged out in the ds1 observation, but there is nothing obviously wrong with that bolometer (it may have much lower high-frequency noise than others).

The simulations actually have more discrepant weights, but that doesn't seem to cause any problems:

The timestreams both have similar artifacts:

while the simulated versions really don't:

This is true even when the relative strength of the atmosphere is higher:

I think the most viable candidate is the 'pointing offset' idea, which will take a little work to simulate properly...

Weighting and Scaling

The "simple" relative sensitivity calibration was causing serious problems.

The assumed model for a "gain" $G$, timestream $S$, and reference timestream $R$ is:
$S = G R$

Naively, one would assume that something like
$G = median(S/R)$
would work. However, it doesn't. The distribution of $G$ values looks like:

Note that this is true no matter what reference bolometer you use . I've monte-carlo tested this and shown that the ratio $S/R<1$ for all bolometers. This is actually a well-understood phenomenon if you look at, e.g., Wall and Jenkins section on linear least squares. If both X and Y values are gaussian-distributed quantities, the linear-least-squares solution is not ideal. Similarly, but more confusing to me, is that the min($\chi^2$) is not a correct solution either, probably because you're minimizing with respect to only one variable.

The solution is Principal Component Analysis. The most correlated component is pretty much the atmosphere (not exactly, but we'll deal with that later). The important thing is that the first row of the eigenvectors gives the relative scaling of each timestream to that most correlated component. Actually the meaning has to be more carefully defined than that, but I can't figure it out right now. However, it looks like the relative scales come out as expected if you scale to a bolometer based on the first row (column?) of eigenvectors corresponding to the most correlated component.

Using PCA and the correct scaling, the results are very pretty, i.e. perfect recovery:

So... now on to the tests of everything ever...

Sampling Causes Problems: Rd 2

Proof that a well-sampled timestream / gaussian is recovered better than a poorly sampled one:

These aren't really the most convincing, since the flux loss is only 1-8% depending on how you count. The atmosphere-included ones have larger flux loss, which is important to understand... WHY does the added noise INCREASE the flux loss?


Ugh, unfortunately, there is a perfect counter-example. The first image shows a timestream in which the input map is sampled by 7.2" pixels, which is just shy of nyquist sampling the 1-sigma width of the gaussian (but is fine for the FWHM of the gaussian). The second shows 3.6" sampling of the same. No improvement. The third, very surprisingly, shows significant improvement - but this was an experiment in which the relative scales were allowed to vary.

The first two each had one high-weight bolometer rejected, the last had 12 high-weight bolos rejected. So that's problably the problem....

Sampling causes problems

One of the simplest experiments that can be run is a point-source observed in point-scan mode, i.e. with smaller step-sizes. I've done this for an Airy disk with noise added in the image plan but not atmosphere (some noise is necessary to avoid bizarre artifacts with weighting when you have truly zero signal and zero noise). At a S/N of 500, the noise is pretty minimal, though.

It turns out for the 'noiseless' images, the PCA cleaning is the issue... curiously, it varies significantly with iteration. Is it worth trying to debug the PCA cleaning for noiseless timestreams? I mean... they really shouldn't have any correlated information anyway.

There is still an outstanding issue where one bolometer gets scaled to be higher than the rest without any apparent reason for doing so.

It is also clear that we do not nyquist-sample the pixels with the timestream. However, that doesn't explain a deficiency observed in the ds1 images. Also note that the sidelobes are not picked up at this S/N, but that's not too surprising.

And here is the explanation: The peak is missed by about 10% because of finite sampling? Somehow the sampled gaussian nearly uniformly underestimates the gaussian... I think this violates theory a bit.....

Here's the problem shown again:

This is a comparison between the input image and a ds1-sampled image with perfectly correlated atmospheric noise. So there is something in the pipeline that is preventing the peaks from achieving the necessary heights.... I wonder if resampling the deconvolved image onto a higher-resolution grid, then downsampling afterwards, would fix this?

Also note that the flux loss increases from 7% to 20% from ds1 to ds5:

Definitive Testing

It is now possible to make artificial timestreams!

First round of tests:
For pointing maps with ~21 scans in 15" steps, array angle is approximately negligible.
deconv does NOT recover point sources, but reconv does (perfectly)
The algorithm starts to decay and produce weird ringey residuals at a S/N ~300. The residuals are only at the 1% level out to S/N ~1000.

Tests running

• Re-doing December 2010 calibration (remember to check & compare AFGL 4029)
  
• Mapping L351 with deconv, reconv, nodeconv
  

ARGH something is wrong with l351 pointing! Is the pointing model broken? What HAPPENED?!


OK, on to 4/8/2011 stuff:
I've done more work on the Memo from a few weeks ago, and now have a concrete recommendation that we do some more observing in 2011. Why I torture myself so, I do not know....

Examination of the 2009 tests shows that reconv is by far the superior approach for pointing observations as previously noted. Since reconv doesn't look bad, but definitely not best, for science, it might be the compromise we have to make. nodeconv completely fails for ds5 data, but the ds5 data in general are abominable. They drop by 15-20% compared to ds1 for Mars. The problem with nodeconv appears to be extremely high noise maps.



...again, forgot to post this on Friday...

Examining deconvolution strategies

This is a direct continuation of what I did yesterday (but only posted today because I forgot to click post).

No-deconvolution doesn't work, at least not reliably. It recovers structures that are too large to be believable. Perhaps a higher threshold should be required to include signal in the model? The noise actually looks too high anyway. Also, the flagging doesn't look great. Grr.

The agreement between the three is somewhat better, down to a 50% increase of reconv over deconv:
[deconv, nodeconv, reconv]:
BMAJ: 0.00916667 BMIN: 0.00916667 PPBEAM: 23.8028 SUM/PPBEAM: 8.19765
Sum: 195.127 Mean: 1.04346 Median: 0.751545 RMS: 0.78342 NPIX: 187
BMAJ: 0.00916667 BMIN: 0.00916667 PPBEAM: 23.8028 SUM/PPBEAM: 10.1592
Sum: 241.816 Mean: 1.29314 Median: 1.04031 RMS: 0.778854 NPIX: 187
BMAJ: 0.00916667 BMIN: 0.00916667 PPBEAM: 23.8028 SUM/PPBEAM: 12.1694
Sum: 289.664 Mean: 1.54901 Median: 1.31107 RMS: 0.808915 NPIX: 187

But this time there is less indication of negative residuals than previously.

I was very confused by negatives being included in the model for nodeconv, then realized it's because of the grow_mask option.

sncut = 1 is now default for nodeconv. I think it makes sense.
(sncut = 1 drops the flux by <10%:
BMAJ: 0.00916667 BMIN: 0.00916667 PPBEAM: 23.8028 SUM/PPBEAM: 9.31507
Sum: 221.725 Mean: 1.18569 Median: 0.915137 RMS: 0.783316 NPIX: 187)

While reconv has produced reasonable results in some cases, a close look at the maps shows that deconv ~ nodeconv << reconv. There is something wrong with reconv. It spreads out and increases the flux artificially. So.... why did it work so damn well for point sources?

The new weighting scheme seems to flag a dangerous number of bolos as 'high weight'. It drops after iteration #1, but not all the way.

Need to remember to reprocess all December 2009 data with more flagged bolos

Deconvolution vs. Not

In the W5 maps, everywhere except AFGL 4029 agrees to within about 20% (better than our supposed offset) independent of deconvolution scheme / modeling scheme, but AFGL 4029 disagrees by a factor of 2-3, indicating a severe dependence on map size.

Here's the AFGL 4029 comparison originally:
BMAJ: 0.00916667 BMIN: 0.00916667 PPBEAM: 23.8028 SUM/PPBEAM: 8.1515
Sum: 194.028 Mean: 1.03758 Median: 0.78441 RMS: 0.760825 NPIX: 187
BMAJ: 0.00916667 BMIN: 0.00916667 PPBEAM: 23.8028 SUM/PPBEAM: 5.77585
Sum: 137.481 Mean: 0.735194 Median: 0.517222 RMS: 0.663245 NPIX: 187
BMAJ: 0.00916667 BMIN: 0.00916667 PPBEAM: 23.8028 SUM/PPBEAM: 14.4043
Sum: 342.863 Mean: 1.83349 Median: 1.62579 RMS: 0.820079 NPIX: 187

Here it is after flagging out another bad bolo:
BMAJ: 0.00916667 BMIN: 0.00916667 PPBEAM: 23.8028 SUM/PPBEAM: 8.14533
Sum: 193.882 Mean: 1.0368 Median: 0.783409 RMS: 0.764047 NPIX: 187
BMAJ: 0.00916667 BMIN: 0.00916667 PPBEAM: 23.8028 SUM/PPBEAM: 10.1518
Sum: 241.642 Mean: 1.2922 Median: 1.06321 RMS: 0.777165 NPIX: 187
BMAJ: 0.00916667 BMIN: 0.00916667 PPBEAM: 23.8028 SUM/PPBEAM: 13.9344
Sum: 331.677 Mean: 1.77367 Median: 1.56277 RMS: 0.820139 NPIX: 187

order is deconvtest, nodeconvtest, 'reconv'

deconvtest has virtually no change. nodeconvtest goes up by a LOT... must be something about the weighting. reconv drops, which is good, but not enough...

curiously, it is not clear that nodeconv converges: [10, 15, 20 iters]
BMAJ: 0.00916667 BMIN: 0.00916667 PPBEAM: 23.8028 SUM/PPBEAM: 8.70939
Sum: 207.308 Mean: 1.34615 Median: 1.09348 RMS: 0.760452 NPIX: 154
BMAJ: 0.00916667 BMIN: 0.00916667 PPBEAM: 23.8028 SUM/PPBEAM: 9.14028
Sum: 217.564 Mean: 1.41275 Median: 1.17223 RMS: 0.759027 NPIX: 154
BMAJ: 0.00916667 BMIN: 0.00916667 PPBEAM: 23.8028 SUM/PPBEAM: 9.42483
Sum: 224.337 Mean: 1.45673 Median: 1.2172 RMS: 0.758978 NPIX: 154

similarly, reconv does not converge:
BMAJ: 0.00916667 BMIN: 0.00916667 PPBEAM: 23.8028 SUM/PPBEAM: 11.5623
Sum: 275.214 Mean: 1.61891 Median: 1.36357 RMS: 0.809771 NPIX: 170
BMAJ: 0.00916667 BMIN: 0.00916667 PPBEAM: 23.8028 SUM/PPBEAM: 12.5848
Sum: 299.553 Mean: 1.76208 Median: 1.51426 RMS: 0.807963 NPIX: 170
BMAJ: 0.00916667 BMIN: 0.00916667 PPBEAM: 23.8028 SUM/PPBEAM: 13.3102
Sum: 316.82 Mean: 1.86365 Median: 1.62092 RMS: 0.804582 NPIX: 170


by contrast, deconv converges rapidly:
BMAJ: 0.00916667 BMIN: 0.00916667 PPBEAM: 23.8028 SUM/PPBEAM: 8.27004
Sum: 196.85 Mean: 0.937381 Median: 0.659373 RMS: 0.759002 NPIX: 210
BMAJ: 0.00916667 BMIN: 0.00916667 PPBEAM: 23.8028 SUM/PPBEAM: 8.39255
Sum: 199.766 Mean: 0.951267 Median: 0.683644 RMS: 0.758766 NPIX: 210
BMAJ: 0.00916667 BMIN: 0.00916667 PPBEAM: 23.8028 SUM/PPBEAM: 8.42881
Sum: 200.629 Mean: 0.955377 Median: 0.690816 RMS: 0.75878 NPIX: 210

while I'm at it, gaussfits for nodeconv and reconv:
Guess: 427,1113 Fit peak: 2.90411 Background: 0.0102637 X,Y position: 426.514822,1113.354863 X,Y FWHM: 64.253564,72.615868 Angle: 0.000000
Guess: 427,1113 Fit peak: 3.03769 Background: 0.000596982 X,Y position: 426.080980,1113.634660 X,Y FWHM: 77.405551,104.292823 Angle: 329.589335
well, isn't that neat! Better than 5% agreement on the peaks. Curiously, the background is higher for nodeconv, which is false: it is directly evident that the background is higher in reconv. So I think this is just a failed fit, unfortunately. The true peak of the reconv AFGL 4029 is 3.9 Jy, so there's a huge residual


Found an error in the noisemap computations that was artificially driving up the noise around NAN pixels in undersampled maps; this probably led to major problems. I fixed it by only "flagging out" pixels with >2 NAN neighbors, i.e. so sparsely sampled that they should be ignored, I hope.... this may have affected modeling in all maps.

Realized the problems with noisemap: if the model exactly equals the weighted mean at a given pixel, the residual will be EXACTLY zero (to within numerical precision). This leads to underestimates of the noise! What we really want is an estimate of the standard deviation on the mean, which is easily computed! Just do a normal weighted sum of the difference between the data and the model (data and the mean); this is also equivalent (conveniently) to a chi^2 statistic, I think... Whoa. How did I not do this before? LETS FIND OUT I'M SURE THERE ARE AWFUL CONSEQUENCES!


Just had another idea to throw in - what if we downweight the scan edges? It won't matter for individual maps or maps of the same size if done uniformly, but it could help incorporate "pointing" maps into "real" maps more accurately.

Deconvolve and Epochs

I've spent a large portion of the last week working on the deconvolver. I found previously that a reconvolved map does a better job of restoring flux than the straight-up deconvolved map for point sources / pointing observations.

However, the same update broke the regular mapping modes, leading to horrible instability in the mapping routines for large maps such as W5. Curiously, it seems that the aspect that breaks is the weighting; somehow the noise drops precipitously in certain bolometers, leading to extremely high weights. Perhaps they somehow dominate the PCA subtraction and therefore have all their noise removed?

Either way, there are a few large-scale changes that need to be made:

1. Since Scaling and Weighting are now done on a whole-timestream basis, we should only map single epochs at once and coadd them after the fact. This approach will also help relieve RAM strain. Since it appears that individual observations are now reasonably convergent with the proper treatment of NANs in the deconvolution scheme, it should be possible to take any individual map and coadd it in a reasonable way.
   
2. Bolometers with bad weights need to be thrown out. Alternatively, and more appropriately, I need to discover WHY their weights are going bad.
   

We also need to explore different weighting schemes.

1. 1/Variance over whole timestream (current default)
   
2. 1/Variance on a per-scan basis (previous default) [based on PSDs]
   
3. Minimum Chi² with Astrophysical Model (??)
   
4. Min Chi² on a per-scan basis?
   

Because of the extensive testing this will require, it is really becoming essential that we develop an arbitrary map creation & testing routine.