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
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:
- 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.
- 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/Variance over whole timestream (current default)
- 1/Variance on a per-scan basis (previous default) [based on PSDs]
- Minimum Chi2 with Astrophysical Model (??)
- Min Chi2 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.
Weird problem
I'm remapping everything, and there's a really strange situation in ic1396... only one field has a source that doesn't have rotation problems, every other observation is clearly improperly rotated. The weird thing is that it's NOT the one you'd expect from the information below: readcol,'/Volumes/disk2/sliced/infiles/ic1396_infile.txt',filelist,format='A'for i=0,6 do begin & ncdf_varget_scale,filelist[i],'array_params',ap & print,filelist[i],ap & endfor/Volumes/disk2/sliced/ic1396/070911_o19_raw_ds5.nc 7.70000 31.2000 70.7000 0.00000 0.00000/Volumes/disk2/sliced/ic1396/070912_o26_raw_ds5.nc 7.70000 31.2000 84.0000 0.00000 0.00000/Volumes/disk2/sliced/ic1396_d/070913_o19_raw_ds5.nc 7.70000 31.2000 84.0000 0.00000 0.00000/Volumes/disk2/sliced/ic1396_l/070913_o17_raw_ds5.nc 7.70000 31.2000 84.0000 0.00000 0.00000/Volumes/disk2/sliced/ic1396_r/070913_o18_raw_ds5.nc 7.70000 31.2000 84.0000 0.00000 0.00000/Volumes/disk2/sliced/ic1396_u/070913_o20_raw_ds5.nc 7.70000 31.2000 84.0000 0.00000 0.00000 One ofthese things is not like the others... but it's 070913_o18_raw_ds5.nc, not /Volumes/disk2/sliced/ic1396/070911_o19_raw_ds5.nc
Progress, but still ds1-ds5 issues
ds1 and ds5 agree pretty well with the recent upgrades to delining and deconvolution. However, there are still counterexamples, e.g. 101208_o13, in which ds5 < ds1:
The 'fitted' values agree better than the 'measured' values now that NANs are treated properly. Spent a few hours today trying to figure out if weighting can explain the difference between ds1 and ds5; it appears to make up for most of it so I'm doing some more experiments. Why is there so much parameter space? Why can't weighting just work? It doesn't.... also wasted a few hours trying to write a python drizzling algorithm, which unfortunately is impossible so I had to resort to an inefficient for loop. Finally got some minor results. It really looks like there is a trend pushing up the recovered flux (i.e. higher volts/Jy) for ds5 over ds1. There is a discrepancy between map types for ds1 but not for ds5, which is actually backwards from what I would have expected, since ds1 will get better-sampled maps.
Luckily, the difference between peak fitting and "measuring" results in very small (insignificant) changes to the calibration curve (recall fitting is direct gaussian fitting; 'measuring' is using the gaussian-fit width and total flux in an ellipse to infer a peak assuming a point source):
Since this work has all been done for the 'bootstrapping' observations that are supposed to tell us if different map sizes are compatible, I have included the map sizes in the diagrams now. However, to really understand the ds1/ds5 difference, there are much better data sets, which I'm now reprocessing using the new and improved methods. (the Whole BGPS is also processing with the new methods in the background, though since the methods are being updated live there may be more changes and it will have to be re-run.... initial looks at W5 are BAD but L030 is GOOD (bordering on amazing))
Careful comparison of ds1 and ds5
I looked very closely at the timestream and maps of 101208_o11 and had a pretty hard time figuring out why the data were different, but it looked like the data really did differ on a point-by-point basis (according to pyflagger). The only conclusion I was able to draw is that the scaling must be off. I realized that the scaling was being done before delining. I moved scaling from readall_pc to premap, and it brought at least this one source into agreement. Time to run ds1-ds5 comparisons again! (this means that ds1 data MUST have deline run on it, but ds5 data doesn't really need it) Here are examples of ds1 and ds5 timestreams, with and without scaling, and ds1 with and without delining:
Trying to bootstrap
I've concluded, based on previous posts http://bolocam.blogspot.com/2011/02/downsampling-why-is-dec-2010-different.html, http://bolocam.blogspot.com/2011/03/revisiting-calibration-yet-again.html, and http://bolocam.blogspot.com/2011/03/workaround-for-individual-maps.html, that ds5 is a problem primarily for undersampled images, i.e. those taken in the normal mapping mode. This makes bootstrapping a bit tricky.
There are two options:
- Map Uranus and AFGL 4029 both in Volts and figure out what flux density AFGL 4029 must have to lie on that curve
- Map Uranus and compute a calibration curve, apply that calibration curve to AFGL 4029, and then compare derived flux densities.
Both have the major problem that the individual AFGL 4029 maps will forcibly be undersampled if I use ds5 data (which is normally OK, according to the first paragraph). In the second case, it is possible to co-add the images and get around the under-sampling issue, while in the first case it is not because of the dependence on total loading (MEANDC).
The real problem is that the whole goal of these observations was to compare the different observing methods and see if they agree (1x1, 3x1, pointing, etc.) since the pointing-style observations were used to calibrate the others. But if the 1x1s are just straight-up unreliable, how can we do the comparison? I think the co-added AFGL 4029 is the only option, but then how do I test if it's correct? It would be really nice to have AFGL 4029 observed with both scan types...
Alright, onto the data. After last week's fix of the bad bolos, I really hope ds1 and ds5 agree. However, first glance at the cal curves says they don't. ds1 and ds2 agree, but ds5 is different.
After checking them out with ds9 *ds[15]*13pca*_map10.fits -scale limits -1 1000 -log -cmap hsv -match colorbars -match scales -match frames wcs &, it appears that the _mask_ data is all... wrong, somehow. That's OK, I want to discard the mask data anyway, so I'm happy to NOT spend time debugging it.
Even after careful examination showing that the fits look good - and noting that the fluxes look pretty much the same - the calibration curves still look rather different. Unfortunately I had to spend 3 hours debugging IDL plotting commands; I want to show the fits each time and save them as postscripts. What does "xyouts" with "/device,/normal" do? I thought that should plot x,y text at the coordinates specified in the plot window... but no, that is JUST /normalize.
Anyway, realized that centroid_map treated NANs as zero. Added ERR keyword (with a reasonable estimate of the error) in centroid_map to ignore NANs. It looks like improper treatment of NANs is responsible for a lot of the scatter seen in the calibration plots.
There is a substantial difference between the "fitted" peak and the "measured" peak (the latter computed by taking the sum of the pixels divided by the area of the fitted gaussian). It looks like the "measured" version is more robust, at first glance. However, unfortunately, for 101208_o11, the difference between ds1 and ds5 exists in both quantities. I will have to examine timestreams now... ARGH.
Well, the timestreams show... that indeed the model is lower in ds1, but not why. The "remainder" (new_astro; the stuff that never gets incorporated into the model but DOES get incorporated into the map) appears to be the same in both. Similarly, there is little to no flux in the PCA atmosphere, so it's not simply being cleaned out. Where is the flux going or coming from?
A workaround for individual maps?
I closely examined the timestreams of 101208_ob7 as I said I would yesterday. Unfortunately, all I can do is describe the symptoms: the first deconvolution model looks good, though it isn't quite as wide as the true source (this should be OK; it is an iterative method, after all). In the second iteration, though, the deconvolution model is even smaller and lower amplitude... and it goes on like that.
Not deconvolving results in a healthy-looking clean map - pretty much what you expect and want to see.
This implies that somehow removing an incomplete deconvolved model leads to more of the source being included in the 'atmosphere' than would have been included with no model subtraction at all. I'm not sure how this is possible. In fact... I'm really quite sure that it is not. The workaround is to only add positive changes to the model. This should 'definitely work' but may be non-convergent and assumes that the model never has anything wrong with it at any iteration. I have demonstrated that this works nicely for the two Uranus observations I tested on, but now I have to run the gamut of tests.... the first (very obvious) problem is that the background is now positive, which is dead wrong. This workaround is not viable. Alright, so what next? I've described the symptoms and that I think they can't occur... A closer look shows that new_astro is not being incorporated into astro_model at the second iteration. Why? AHA! Pyflagger + find_all_points reveals the problem!
Map value: 16.939728 Weighted average: 17.476323 Unweighted Average: 524.573136 scan,bolo,time: mapped astro flags weight scale 3, 22, 12: 8.380408 13.561113 0.000000 0.025132 1.000000 4, 124, 23: 822.005327 13.561113 0.000000 0.000038 1.118012 4, 21, 38: 719.408983 13.561113 0.000000 0.000037 0.946721 5, 20, 7: 4.470616 13.561113 0.000000 0.013303 1.400000 5, 119, 23: 882.508303 13.561113 0.000000 0.000033 0.926887 5, 100, 35: 327.007750 13.561113 0.000000 0.000074 1.184397 5, 106, 38: 162.562098 13.561113 0.000000 0.000704 0.970000 6, 116, 27: 779.075640 13.561113 0.000000 0.000033 0.891768 8, 112, 3: 235.557390 13.561113 0.000000 0.000147 0.947130 9, 3, 14: 966.721773 13.561113 0.000000 0.000032 1.166292 9, 109, 41: 139.753656 13.561113 0.000000 0.000753 1.075269 10, 104, 8: 641.121935 13.561113 0.000000 0.000050 0.927827 10, 105, 24: 4.323228 13.561113 0.000000 0.032759 0.019022 10, 32, 36: 847.646990 13.561113 0.000000 0.000034 1.099406 11, 36, 9: 834.757586 13.561113 0.000000 0.000038 1.184751 11, 76, 37: 566.851891 13.561113 0.000000 0.000040 1.111000 12, 77, 13: 834.603090 13.561113 0.000000 0.000034 1.128464 12, 44, 44: 335.465654 13.561113 0.000000 0.000195 2.165775 13, 26, 17: 50.423143 13.561113 0.000000 0.004826 0.829932 13, 75, 29: 724.884676 13.561113 0.000000 0.000042 0.923077 14, 49, 21: 797.618990 13.561113 0.000000 0.000038 1.091918 14, 29, 33: 743.856012 13.561113 0.000000 0.000035 1.050360 15, 33, 13: 660.670099 13.561113 0.000000 0.000031 0.832180 15, 53, 25: 604.174286 13.561113 0.000000 0.000047 0.889922 15, 88, 40: 4.626476 13.561113 0.000000 0.008241 0.191489 17, 64, 20: 778.950533 13.561113 0.000000 0.000037 1.233108 18, 68, 30: 686.048136 13.561113 0.000000 0.000040 1.387283
Note that the lowest points have the highest weights. They DEFINITELY shouldn't. What's wrong with them? Apparently they have NO sensitivity to the sky! What?! There were a bunch of bad bolos in Dec2010 that weren't flagged out... I wonder if that problem persists to other epochs. Still, why does it only affect pointing observations? Looking at the power spectra... the large-timescale stuff becomes less dominant when scans are longer, but the noisy spectra are still clearly noise-only. How odd. Dropped to 112 good bolos from 134. That is much more believable. Have to go back and fix Dec09 data too... Even after fixing the bad bolos, the model drops with iteration number. Why why why? Well, looking at deconv_map, I've always returned the truly deconvolved version, not the reconvolved... maybe the reconvolved really is better? Again, this will have to be extensively tested, but it certainly gets rid of the obvious/dominant error that the model kept dropping off. However, FINALLY, based on how ridiculously good the reconv-deconvolved map looks, I think I'm ready to do the extensive pipeline tests. So, 10dec_caltest has been started up with all of the new bolo_params applied and the changes in place to deconv_map... let's see what happens.
After that runs, I'll have to re-run the fit_and_plot routines
Revisiting calibration yet again
The recent hiatus for paper revisions has, unfortunately, come to an end. Re-examining my work, I did quite a lot but encountered many dead-ends. First, we would very much like to use an *identical* reduction process on both the calibration data and the science data. That way, we could feel very confident that the reduction process isn't introducing any weird artifacts. Unfortunately, I discovered early on that using ds5 data, 13 pca components, and n>1 iterations resulted in strange shape and flux conservation failures. These errors do NOT occur in co-added maps; they are unique to single-observation scans (though I don't recollect whether 2 scans is enough or if you need more). I spent many hours banging my head against this problem and have never gotten a satisfactory solution. But perhaps it's time to approach it again. The map00 images look MUCH rounder and generally better than the map10 images. So, the problem I need to examine is the iterative process. Why does it fail for single images? Is it something about the noise properties? model00 looks fine... what gets put into the timestream? Examining timestreams is a terrible and horrendous process... but what else can I do? The next step will be to examine the timestreams of a particular observation. I think a good choice is 101208_ob7; the next observation, 101208_ob8 was a large-area map and it looks fine (i.e., it improves with iteration). So I can start looking at the effects of polysub, iteration, etc. on this particular source. Of course, the stupid trick with the pipeline - every time - is that "fixing" a problem for one source has a nasty tendency to break it for all other sources. That's why there are so many flags that can be passed around. Still, this is the approach I have to take...
Downsampling - Why is Dec 2010 different?
As you'll recall from my previous post (and references therein...), the 2005 Orion data shows a discrepancy between ds1 and ds5 data in which the ds1 data is significantly (~10%) higher than the ds5 data. However, the 2010 Uranus observations show much larger discrepancies between ds1 and ds5 favoring the ds5 data! Because that was somewhat unbelievable to me, I ran ds1-ds5 comparisons on Uranus data from other epochs, and discovered that ds1>ds5 uniformly (also, it looks a LOT better). So, the question remains, WHY is the Dec 2010 data brighter in ds5? More confusing to me, why do the ds5 PSFs from 2010 look so reasonable, while the ds5 PSFs from all earlier epochs look terrible? For example, I use the observations 070727_o31 and _o32. These show clearly the blurring and flux-loss that happens when ds5 data are used for 'normal' point-sources:
Downsampling - what is going on?
The downsampling failure I noted previously appears to be illusory. It may be that the offset noted only holds for single-frame images, in which there may be many blank pixels. It is possible - though not certain - that the ds1 images were significantly higher than ds5 because more noise-only pixels were included with higher outliers; i.e., ds1 high-outlier noise was being compared to ds5 noise that was lower amplitude.
What led to these conclusions? First, I was getting inconsistent results looking at Uranus in particular - ds5 appeared to have higher fluxes than ds1. This was inconsistent with earlier results on OMC1. Partly, this is because I switched from my hacked-together plots to the much more refined compare_images script, which demonstrated the effect of changing the cutoff of the comparison.
Also, I added in a Pearson Correlation Coefficient computation. Given a single data set with the only difference being downsampling, the data should be perfectly correlated even if there is a flux offset (correlation should be 1, but the best fit slope should not be). It was an indication of a problem when I started seeing correlation coefficients <0.90 for data that had already been sigma-cut; that means that noise was being included in the correlation computations.
Therefore, the approach needed is to cut out the high pixels that are on map edges. This I accomplished by adding an 'aperture' capability to the compare_images code (for Uranus) and cropping using montage and a wcs-based box for Orion.
The results... are ambiguous. Wow. In some sub-fields - within the same co-added map - the agreement is near-perfect.
In others, ds1 is clearly > ds5.
What's going on? ds1 does look uniformly more smooth. Note that the disagreement is nearly scale-free:
OK, so given the conclusion in Orion that ds1>=ds5, what's the deal with Uranus?
The first two comparisons are for 1x1° observations; in both cases ds1 < ds5, but by 6% and 24% respectively! The image of Uranus looks much better (because of lack of parallel lines) in the second, more extreme case. In both cases, the ds5 excess is nearly scale-free (not shown).
The 3x1s are also highly discrepant. #12 shows nearly perfect agreement, albeit with high dispersion (low correlation) because of pixel-to-pixel variations around the peak. #13 is the only observation with a huge DS1 excess. It also demonstrates very poor correlation. It looks like the telescope got bumped for the ds5 data (which is not actually possible; recall they're the same data set). What happened here? Maybe a glitch that went unflagged (mad_flagger is off by default for individual scans)?
In observations 4 and 5, we're looking at a 40-50% excess in ds5! What the heck? There really is no clear explanation for this.
But... what? Magically, they come into perfect agreement when the scan axis nearly lines up with the coordinate axis! Or, is this just an effect of the worse weather on night 2?
Next thing to try: masked source map comparison. Unfortunately, masking royally screwed up the long scans - probably because the initial polysub didn't work. And masking in the individual point source maps did nothing... so that pretty much rules out atmospheric oversubtraction, doesn't it? What else could be causing this offset? 0pca looks the same as 13pca, give or take, so it's not the atmospheric subtraction. Could the downsampling result in an offset in the bolo-scaling? Where else in the process could things go wrong? Tomorrow, need to investigate .sav files with pyflagger...
