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

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.

Weighting and high-frequency noise

Image of PSDs (with no normalization) of the raw (blue), delined and exponential and polynomial subtracted (white), the noise timestream (yellow), and the data (cyan).

The good: It looks like all of the powerline noise got into the noise timestream and almost none in the data.

The bad: weighting is based on the noise timestream so it's possible that the weights aren't quite right as a result

The weird: the data PSD. What's up with that? Apparently I'm preferentially subtracting certain scales but I don't know why, unless deconvolution is at fault.


Edit/Update: The deconvolution is definitely at fault. Here's the same scan done without deconvolution:



It should have been obvious; the cyan in the first plot is the PSD of the deconvolution straight up, and that should have no high-frequency structure...

Implementing weighting

Not as easy as it ought to be.


I think I need to do a few things:
1. check and make sure there are no more of those !@#$!@#$#@% different sized array subtractions/multiplications. 'weight' and 'best_astro_model' need to have the same size & shape in mem_iter_pc
2. I guess just check and make sure stuff works. The weighted mean I'm using appears to be right: sum(weight * value) / sum(weight)

I hate making lists that end up being two items....