Deconvolution vs. Not

I've done some by-eye comparisons of deconvolutions vs no deconvolution. The no-deconvolution clearly does better on the bright sources: probably the deconvolved beam size is a little bit different from (larger? smaller?) the actual source size. Deconvolution does better in putting noise from noisy regions into the noisemap and keeping it out of the astromap.

Both might be useful - the deconvolved maps may end up being prettier, but the no-deconvolve maps will probably have more reliable fluxes.

This all probably relies on testing / simulation.

4 errors

l047, l083, l359 (l000 is in /usb, l359 is in /scratch), l136p15

Only 83 is a code error, the others are file location errors. l083 is an error in 'pixshift' with subscripting and will require real debugging.

def_user_common modified

The NP doesn't work without USER_COMMON defined because at a single point it makes use of the 'which' function.

def_user_common is a disaster, though, because it calls get_screen_size(), which forces an X connection to be launched, which means that if you start a remote session and close X IDL crashes. This has been a constant nagging problem and has cost me ~a few days of work. So I commented out the offending (and offensive) line. The worst of it is, I'm not even convinced there's anything that uses any of that information. The common block is only needed to get paths, so which could have been written better.

The pipeline ground to a halt for an unknown reason at an unknown point so I added a lot more timing outputs to try to figure out what's going on.

Back to the grind...

flagger

Discovered that the flagger causes some pretty serious problems if you try to run it on a non-coadd for obvious reasons: there's practically no coverage! Individual maps should NOT go through the flagger, so I have added a piece of code that turns off the flagger if only one observation is being used. That code is in map_iter where the flagger is called.

Version 0.6 begun

I started on (what I hope will be) v0.6 tonight.

New things:
-Sigma-rejection flagging in the spatial domain
-Inverse-variance weighting

MAD flagger

; The Mad Flagger
; Flag based on the median average deviation within a spatial pixel
function mad_flagger,data,inds,flags,nsig=nsig

    t = systime(/sec)
    f0 = total(where(flags))

    if n_e(nsig) eq 0 then nsig = 3    
    
    newflags = flags

    mx=max(inds)
    vec3=fltarr(mx+1)
    h=histogram(inds,reverse_indices=ri,OMIN=om)
    for j=0L,n_elements(h)-1 do begin
        if ri[j+1] gt ri[j] then begin
            v_inds = [ri[ri[j]:ri[j+1]-1]]
            if n_e(v_inds) gt 2 then begin
                vec = data[v_inds]
;                vecmad = mad(vec)  ; the MAD is WAY too small!  I ended up rejecting 8% of points!
                vecmad = stddev(vec)
                vecmed = median(vec,/even)
                madreject = where( (vec gt vecmed + nsig*vecmad) or (vec lt vecmed - nsig*vecmad) )
                if (n_e(madreject) gt 0 and total(madreject)) gt 0 then begin
                    reject_inds = v_inds[madreject]
                    newflags[reject_inds] = 1
                endif 
            endif
        endif
    endfor
    print,"MAD flagger took ",strc(systime(/sec)-t)," seconds and flagged ",$
        strc(round(total(where(newflags)) - f0)),' points'

    return,newflags

end

Median drizzling

http://groups.google.com/group/comp.lang.idl-pvwave/browse_thread/thread/762770933591238d

Someone going by the name 'wox' introduced me to a spectacularly simple drizzling algorithm for medianing. It will be a little less efficient because it can't take advantage of some 'partial sum' tricks that the average can, but it will median over ALL data points, which is an advantage.

Massive upgrades, new strategy?

A few major things today:
1. I renamed 'ts_to_map' to 'drizzle' and made 'ts_to_map' a wrapper that will either median combine by scans (if passed the right keywords, i.e. tstomapmedian and scans_info). This will reject cosmic rays much more effectively than previously.
2. My extensive test runs on the 'weird' fields completed. Conclusions:
a. Deconvolution is DEFINITELY responsibly for the weird fuzzy effects seen esp. in l=33, l=2.
b. Even a descending iterator is ineffective at alleviating this problem
What I still don't understand is WHY the deconvolver is behaving badly for some regions. It is extracting peaks that are too high.
c. Not deconvolving works pretty well, so I'm not really concerned. I'm resetting the defaults to no deconvolution once my next test is done....
3. Started up a new test run to compare median stacking with weighted average stacking

5x1?

060609_o17 is weird looking. It appears to be TWO 3x1s with an overlap, but it's only 1 observation. Any idea what happened?

Weighting - next steps

Idea: 'edge taper' for scans. The scan start/end usually has higher noise, so it would be good to downweight those regions slightly. Probably only the first/last 5% of each scan should be linearly downweighted.

I'm running a lot of tests on l001 with/without deconvolution, lots of iterations, different pca components. They'll be done sometime tomorrow.

Paragraph for the methods paper

The raw data from Bolocam contains noise components from the atmosphere and instrument in addition to the astrophysical signal. To remove the atmospheric noise, an iterative approach was required.

1. The median across all bolometers is subtracted
   
2. A set number of principle components are subtracted. The principle components are the most correlated components between bolometers. In this process both the atmosphere above the telescope - which is assumed to be constant across the field of view - and any large-scale astrophysical structure are removed.
   
3. The timestream data is mapped into the plane of the sky. Data points are mapped to the nearest pixel. 7.2" pixels are used so that sampling is better than Nyquist.
   
4. The map is deconvolved using a maximum entropy deconvolution algorithm ( Based on paper by Hollis, Dorband, Yusef-Zadeh, Ap.J. Feb.1992, written by Frank Varosi at NASA/GSFC 1992)
   
5. The deconvolved map is returned to a timestream and subtracted from the original to yield a noise-only timestream.
   
6. Power spectral densities are calculated for each scan in the noise timestream, and weights are calculated from these. [At the moment, the weights are actually inverse-variance]
   
7. The deconvolved map timestream is subtracted from the raw timestream, and then steps 1-6 are repeated on that timestream to recover flux that was oversubtracted in the first iteration.
   

Convergence takes ???? iterations....
??? PCA components are subtracted [default 13]...

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

Rewriting pipeline again

Not entirely, but I realized my code had become... clunky, at best.

I'm rewriting the wrappers using structs. I HATE IDL structs, but they're necessary.... if only they would make them dynamically modifiable.

Weighting worked to some degree, the mapping is pretty much done.

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

Latest re-run of the pipeline

Wow. A lot of errors this time around.


bash: /scratch/adam_work/l020/l020_fitslist.txt: Permission denied
MRDFITS: Image array (1634,560) Type=Real*4
% HEULER: FITS header is already in Galactic: nothing changed
% READFITS: ERROR - Unable to locate file
060628_o32_raw_ds5.nc_indiv13pca_map01.fits
% Variable is undefined: HDIMAGE.
% Execution halted at: IMAGE_SHIFTS 44
/home/milkyway/student/ginsbura/bgps_pipeline/alignment/image_shifts.pro
% COALIGN_FIELD 31
/home/milkyway/student/ginsbura/bgps_pipeline/wrappers/coalign_field.pro
% $MAIN$


Updated radec_offsets in /scratch/sliced_polychrome/l034/050708_o21_raw_ds5.nc
Writing offsets to -0.00192975 -9.22609e-09 (ncfile was -0.000337092)
% ALIGNMENT_OFFSETS_TO_NCDF: Could not find file
% Execution halted at: ALIGNMENT_OFFSETS_TO_NCDF 4
/home/milkyway/student/ginsbura/bgps_pipeline/preproc/alignment_offsets_to_
ncdf.pro
% WRITE_IMSHIFTS 11
/home/milkyway/student/ginsbura/bgps_pipeline/alignment/write_imshifts.pro
% COALIGN_FIELD 34
/home/milkyway/student/ginsbura/bgps_pipeline/wrappers/coalign_field.pro
% IMAGE_SHIFTS 44
/home/milkyway/student/ginsbura/bgps_pipeline/alignment/image_shifts.pro
% COALIGN_FIELD 31
/home/milkyway/student/ginsbura/bgps_pipeline/wrappers/coalign_field.pro
% $MAIN$


/scratch/adam_work/l354/070709_o16_raw_ds5.nc_indiv13pca took 110.42429 sec.

MRDFITS: File access error
% HEULER: ERROR - First parameter must be a FITS header or astrometry
structure
% Execution halted at: HEULER 78

FIELD l351 COMPLETED at Sat Nov 8 05:34:01 2008. Took 1122.3564 seconds
% Program caused arithmetic error: Floating underflow
% Program caused arithmetic error: Floating illegal operand
FIELD l052 BEGUN at Sat Nov 8 05:34:01 2008
% FILE_LINES: Error opening file.
File: /scratch/sliced/l052/l052_infile.txt
No such file or directory
% Error occurred at: READCOL 141

FIELD l359 COMPLETED at Sat Nov 8 05:41:34 2008. Took 452.60051 seconds
% Program caused arithmetic error: Floating underflow
% Program caused arithmetic error: Floating illegal operand
FIELD l066 BEGUN at Sat Nov 8 05:41:34 2008
% FILE_LINES: Error opening file.
File: /scratch/sliced/l066/l066_infile.txt
No such file or directory
% Error occurred at: READCOL 141


FIELD l072 BEGUN at Sat Nov 8 05:41:34 2008
% FILE_LINES: Error opening file.
File: /scratch/sliced/l072/l072_infile.txt
No such file or directory
% Error occurred at: READCOL 141



MRDFITS: Image array (628,561) Type=Real*4
% HEULER: FITS header is already in Galactic: nothing changed
% HEULER: FITS header is already in Galactic: nothing changed
/scratch/adam_work/l086/060627_o23_raw_ds5.nc_indiv13pca_map01.fits
shifts in pixels: -9.77712 11.8034 shifts in arcsec l/b:
70.395277 84.984804 delta-ra (as): 19.007974
delta-dec (as): -108.70337 x/y error in pixels: 3.99631
0.0701244 ra/dec error: -0.0079926204 0.00014024887 gal l/b error:
-0.0079926204 0.00014024887
% HEULER: FITS header is already in Galactic: nothing changed
WARNING! x,y offs were too large: 203.757 -0.224758
so a different method is being attempted
WARNING! x,y offs were too large on the SECOND TRY too: -74.7339
-159.448 so a third method is being attempted
This third method is meant to deal with small overlap regions that may be confused by high noise, it restricts
to a maximum 15 pixel offset (really less than that in practice) and therefore should be checked carefully.
/scratch/adam_work/l086/060627_o26_raw_ds5.nc_indiv13pca_map01.fits
shifts in pixels: -1.38579 1.39157 shifts in arcsec l/b:
9.9777001 10.019338 delta-ra (as): 1.1521393
delta-dec (as): -14.093081 x/y error in pixels: 0.0814682
0.0882818 ra/dec error: -0.00016293635 0.00017656364 gal l/b error:
-0.00016293635 0.00017656364
% HEULER: FITS header is already in Galactic: nothing changed
WARNING! x,y offs were too large: 9.13210 11.3905
so a different method is being attempted
WARNING! x,y offs were too large on the SECOND TRY too: -79.0260
20.0907 so a third method is being attempted
This third method is meant to deal with small overlap regions that may be confused by high noise, it restricts
to a maximum 15 pixel offset (really less than that in practice) and therefore should be checked carefully.
/scratch/adam_work/l086/070707_o49_raw_ds5.nc_indiv13pca_map01.fits
shifts in pixels: -9.13210 -11.3905 shifts in arcsec l/b:
65.751123 -82.011258 delta-ra (as): -105.06683
delta-dec (as): 3.1919479 x/y error in pixels: 0.486307
0.950778 ra/dec error: -0.00097261386 0.0019015568 gal l/b error:
-0.00097261386 0.0019015568
% HEULER: FITS header is already in Galactic: nothing changed
WARNING! x,y offs were too large: 0.384630 13.1890
so a different method is being attempted
WARNING! x,y offs were too large on the SECOND TRY too: 205.434
-23.9600 so a third method is being attempted
This third method is meant to deal with small overlap regions that may be confused by high noise, it restricts
to a maximum 15 pixel offset (really less than that in practice) and therefore should be checked carefully.
/scratch/adam_work/l086/070708_o37_raw_ds5.nc_indiv13pca_map01.fits
shifts in pixels: -0.384630 -13.1890 shifts in arcsec l/b:
2.7693373 -94.960636 delta-ra (as): -74.063617
delta-dec (as): 59.509407 x/y error in pixels: 2.29135
0.933383 ra/dec error: -0.0045826990 0.0018667649 gal l/b error:
-0.0045826990 0.0018667649
% HEULER: FITS header is already in Galactic: nothing changed
WARNING! x,y offs were too large: -143.139 10.8130
so a different method is being attempted
% Subscript range values of the form low:high must be >= 0, < size, with low
<= high: CORR.
% Execution halted at: PIXSHIFT 107




ALL PREPROC ... output is /scratch/adam_work/l133/070911_o25_raw_ds5.nc_indiv13pca READING IN 1 FILES FROM /scratch/sliced/l133/070911_o25_raw_ds5.nc
Reading file /scratch/sliced/l133/070911_o25_raw_ds5.nc
Applying pointing model with alt/az off -0.024687220 -0.029154389
NO OFFSETS APPLIED
3.6277909 sec.
No values flagged as bad!
AC_BOLOS has 0.0 NAN points, which is 0.0 of total
Flatfielding (polynomial subtraction) ...0.92621708 sec.
No values flagged as bad!
AC_BOLOS has 0.0 NAN points, which is 0.0 of total
Prepare Map ... 1.4727209 sec.
% CHECK_FITS: NAXIS keywords in FITS header have been updated
Done preproc ... 11.113711 sec.
ITERATION NUMBER 0 ...SKY SUBTRACTION - MEDIAN0.15759706 sec.
BEGINNING PCA SUBTRACTION OF COMPONENTS 0 TO 13 ... 0.68130398 sec.
Flatfielding (polynomial subtraction) ...0.96450186 sec.
% CHECK_FITS: NAXIS keywords in FITS header have been updated
% CHECK_FITS: NAXIS keywords in FITS header have been updated
% Program caused arithmetic error: Floating underflow
Deconvolution iteration 0: 1.2357290 sec.
Deconvolution iteration 1: 0.93198013 sec.
Deconvolution iteration 2: 0.93152189 sec.
Deconvolution iteration 3: 0.93164992 sec.
Chi2 for iteration 0 is 1.22593e+07 with 2.07110e+06 degrees of freedom
FINISHED ITERATION 0 13.803532 sec.
ITERATION NUMBER 1 ...SKY SUBTRACTION - MEDIAN0.15691710 sec.
BEGINNING PCA SUBTRACTION OF COMPONENTS 0 TO 13 ... % TRIQL: Too many iterations in tqli.
% Error occurred at: EIGENQL 122
/usr/local/itt/idl70/lib/eigenql.pro
% PCA_SUBTRACT 19



/scratch/adam_work/l135p1/060917_o12_raw_ds5.nc_indiv13pca took 30.478553 sec.

MRDFITS: File access error
% HEULER: ERROR - First parameter must be a FITS header or astrometry
structure
% Execution halted at: HEULER 78


MRDFITS: Image array (1630,555) Type=Real*4
% HEULER: FITS header is already in Galactic: nothing changed
% HEULER: FITS header is already in Galactic: nothing changed
/scratch/adam_work/l068/070910_ob3_raw_ds5.nc_indiv13pca_map01.fits
shifts in pixels: -0.0271645 -0.0270424 shifts in arcsec l/b:
0.19558409 -0.19470519 delta-ra (as): -0.26866187
delta-dec (as): -0.062403024 x/y error in pixels: 0.000193238
0.00126581 ra/dec error: -3.8647498e-07 2.5316247e-06 gal l/b error:
-3.8647498e-07 2.5316247e-06
% HEULER: FITS header is already in Galactic: nothing changed
WARNING! x,y offs were too large: 34.8510 -5.72119
so a different method is being attempted
WARNING! x,y offs were too large on the SECOND TRY too: 142.704
76.4267 so a third method is being attempted
This third method is meant to deal with small overlap regions that may be confused by high noise, it restricts
to a maximum 15 pixel offset (really less than that in practice) and therefore should be checked carefully.
/scratch/adam_work/l068/070911_ob8_raw_ds5.nc_indiv13pca_map01.fits
shifts in pixels: 15.0101 3.42143 shifts in arcsec l/b:
-108.07276 24.634271 delta-ra (as): 78.221561
delta-dec (as): 78.545365 x/y error in pixels: 0.0114306
0.184410 ra/dec error: -2.2861189e-05 0.00036881918 gal l/b error:
-2.2861189e-05 0.00036881918
% HEULER: FITS header is already in Galactic: nothing changed
WARNING! x,y offs were too large: 56.2296 14.6994
so a different method is being attempted
WARNING! x,y offs were too large on the SECOND TRY too: 728.753
-184.087 so a third method is being attempted
This third method is meant to deal with small overlap regions that may be confused by high noise, it restricts
to a maximum 15 pixel offset (really less than that in practice) and therefore should be checked carefully.
/scratch/adam_work/l068/070911_ob9_raw_ds5.nc_indiv13pca_map01.fits
shifts in pixels: -56.2296 -14.6994 shifts in arcsec l/b:
404.85292 -105.83592 delta-ra (as): -304.51436
delta-dec (as): -286.88182 x/y error in pixels: 7.32951
0.00354093 ra/dec error: -0.014659024 7.0818565e-06 gal l/b error:
-0.014659024 7.0818565e-06
% HEULER: FITS header is already in Galactic: nothing changed
WARNING! x,y offs were too large: 3.31224 -13.8321
so a different method is being attempted
% Subscript range values of the form low:high must be >= 0, < size, with low
<= high: CORR.
% Execution halted at: PIXSHIFT 107

/scratch/adam_work/l111/060605_o48_raw_ds5.nc_indiv13pca took 10.098406 sec.

MRDFITS: File access error
% HEULER: ERROR - First parameter must be a FITS header or astrometry
structure
% Execution halted at: HEULER 78


Writing offsets to /scratch/sliced/l027/060616_o10_raw_ds5.nc 0.00120866
-0.00287235 (ncfile was 060616_o10_raw_ds5.nc)
Updated radec_offsets in /scratch/sliced/l027/060616_o10_raw_ds5.nc
Writing offsets to 5.77312e-05 -2.39760e-06 (ncfile was 0.000409983)
% ALIGNMENT_OFFSETS_TO_NCDF: Could not find file
% Execution halted at: ALIGNMENT_OFFSETS_TO_NCDF 4


BGPS Photo

http://www.nrao.edu/index.php/learn/gallery/imagecontest