Adam Ginsburg's blog

• 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.

---

• 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:

  1. Map Uranus and AFGL 4029 both in Volts and figure out what flux density AFGL 4029 must have to lie on that curve
  2. 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 - 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...

---

• Making maps faster

  The fundamental problem at this point is making the pipeline run faster. At current speeds, with undownsampled data, it may take ~days to process a single map. Ideas for faster processing:

  1. Find out how long it takes to converge to 1%, 5%.... If it only requires 10 iterations, that's a factor of 2 savings over current strategies.
  2. Use downsampled data of some sort if possible. Does DS2 match DS1? How do we measure it? Flux-flux comparison and PSF point-source size measurements are the most important. Need to automate PSF comparison....
  3. Can we use median-combined individual images as a 0th order model? I bet the answer is 'yes' and will probably increase the speed of convergence by a large amount. Tests to run? This is probably needed if we are to split up the 'super-fields' into smaller sub-fields, otherwise overlapping data will be used less effectively.
  4. Find some way to keep bgps.raw, bgps.ra, bgps.dec, and other items that are only used once on the HD during the iterative process. Is there any way to separate out data in a struct in this manner?

---

• Delining - Maps

  First comment - delining has no effect on downsampled data. At least for the 0709 epoch, there were NO lines AT ALL in the data. From 0-5 Hz, it was just empty. So we don't have to worry about that... the problem only affects fully-sampled data. Then, onto map comparisons. Curiously, the noise levels don't drop after delining. They actually go up a bit. This may be because of the effects on PCA cleaning. However, flux levels in the sources go up by 0-10%. As usual, the change in flux changes from field to field without any obvious reason. Example 1: A pointing field. The source is ~2% brighter in the delined version, but otherwise the match between the two is nearly perfect.

  

  Example 2: A bigger map, where the flux recovery is much greater when delining, but the background levels are also higher.

  

---

• Delining and the Cleaning process

  One item I forgot to mention last night was the effects of lines/delining on PCA subtraction. These should be the primary effects on the final map for all epochs except 2010, in which case the primary effect SHOULD be to reduce substantial noise. In the examples below, there are PSDs of whole timestreams (left) and example timestreams from single scans (right). The first thing to note is that the delined timestreams still have correlated components in the line region, but they are suppressed - their amplitudes, and therefore their sort order in the PCA removal scheme, are changed. Since PCA cleaning is by its nature adaptive (the number of components remains fixed, but the order changes), these effects can be significant and dangerous. If the line noise is more correlated, a PCA component will be dedicated to removing it instead of atmospheric signal. Below are examples from l089 (epoch 0709) first. These have less correlated line noise and are more typical of BGPS observations. The first PCA component, the average, does not change much with PCA cleaning. However it is clear that the second component changes substantially, from large-amplitude high-frequency noise to small-amplitude variations that are very likely to describe atmosphere.

  ---

  Example 1 - Zeroing out the lines:

  

  Example 2 - Fitting and removing the lines:

  

  Example 3 - Suppressing the lines with a non-fitted Gaussian:

  

  ---

  The next examples are from December 2010 observations of Uranus. In this case, the correlated noise component is clearly dominant.

  Zeroing lines:

  

  Fitted lines:

  

  Non-fitted gaussian suppression:

  

  ---

  Finally, these two are demonstrations of what you might expect to see for a purely noiseless images of a planet (it was constructed from a PSF). PCA is first:

  

  A single bolometer's timestream and PSD:

  

---

« Page 2 / 5 »