forked from ACTCollaboration/enlib
-
Notifications
You must be signed in to change notification settings - Fork 0
/
todfilter.py
136 lines (129 loc) · 5.54 KB
/
todfilter.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
import numpy as np, time, h5py
from enlib import config, fft, utils, gapfill, todops, pmat, rangelist
config.default("gfilter_jon_naz", 8, "The number of azimuth modes to fit/subtract in Jon's polynomial ground filter.")
config.default("gfilter_jon_nt", 10, "The number of time modes to fit/subtract in Jon's polynomial ground filter.")
config.default("gfilter_jon_nhwp", 0, "The number of hwp modes to fit/subtract in Jon's polynomial ground filter.")
config.default("gfilter_jon_niter", 3, "The number of time modes to fit/subtract in Jon's polynomial ground filter.")
def filter_poly_jon(tod, az, weights=None, naz=None, nt=None, niter=None, cuts=None, hwp=None, nhwp=None, deslope=True, inplace=True):
"""Fix naz Legendre polynomials in az and nt other polynomials
in t jointly. Then subtract the best fit from the data.
The subtraction is inplace, so tod is modified. If naz or nt are
negative, they are fit for, but not subtracted.
NOTE: This function may leave tod nonperiodic.
"""
#moomoo = tod[:8].copy()
naz = config.get("gfilter_jon_naz", naz)
nt = config.get("gfilter_jon_nt", nt)
nhwp= config.get("gfilter_jon_nhwp", nhwp)
niter = config.get("gfilter_jon_niter", niter)
if not inplace: tod = tod.copy()
do_gapfill = cuts is not None
#print "Mos", naz, nt, nhwp
#print hwp
# No point in iterating if we aren't gapfilling
if not do_gapfill: niter = 1
if hwp is None or np.all(hwp==0): nhwp = 0
naz, asign = np.abs(naz), np.sign(naz)
nt, tsign = np.abs(nt), np.sign(nt)
nhwp,hsign = np.abs(nhwp),np.sign(nhwp)
d = tod.reshape(-1,tod.shape[-1])
if naz == 0 and nt == 0 and nhwp == 0: return tod
# Build our set of basis functions. These are shared
# across iterations.
B = np.zeros([naz+nt+nhwp,d.shape[-1]],dtype=tod.dtype)
if naz > 0:
# Build azimuth basis as polynomials
x = utils.rescale(az,[-1,1])
for i in range(naz): B[i] = x**(i+1)
if nt > 0:
x = np.linspace(-1,1,d.shape[-1],endpoint=False)
for i in range(nt): B[naz+i] = x**i
if nhwp > 0:
# Use sin and cos to avoid discontinuities
c = np.cos(hwp)
s = np.sin(hwp)
for i in range(nhwp):
j = i/2+1
x = np.cos(j*hwp) if i%2 == 0 else np.sin(j*hwp)
B[naz+nt+i] = x
for it in range(niter):
if do_gapfill: gapfill.gapfill(d, cuts, inplace=True)
# Solve for the best fit for each detector, [nbasis,ndet]
# B[b,n], d[d,n], amps[b,d]
if weights is None:
try:
amps = np.linalg.solve(B.dot(B.T),B.dot(d.T))
except np.linalg.LinAlgError as e:
print "LinAlgError in todfilter. Skipping"
continue
else:
w = weights.reshape(-1,weights.shape[-1])
amps = np.zeros([naz+nt+nhwp,d.shape[0]],dtype=tod.dtype)
for di in range(len(tod)):
try:
amps[:,di] = np.linalg.solve((B*w[di]).dot(B.T),B.dot(w[di]*d[di]))
except np.linalg.LinAlgError as e:
print "LinAlgError in todfilter di %d. Skipping" % di
continue
#print "amps", amps[:,2]
# Subtract the best fit
if asign > 0: d -= amps[:naz].T.dot(B[:naz])
if tsign > 0: d -= amps[naz:naz+nt].T.dot(B[naz:naz+nt])
if hsign > 0: d -= amps[naz+nt:naz+nt+nhwp].T.dot(B[naz+nt:naz+nt+nhwp])
# Why was this necessary?
if do_gapfill: gapfill.gapfill(d, cuts, inplace=True)
if deslope: utils.deslope(tod, w=8, inplace=True)
res = d.reshape(tod.shape)
return res
def deproject_vecs(tods, dark, nmode=50, cuts=None, deslope=True, inplace=True):
"""Given a tod[ndet,nsamp] and a set of basis modes dark[nmode,nsamp], fit
each tod in to the basis modes and subtract them from the tod. The fit
ignores the lowest nmode fourier modes, and cut regions are approximately ignored."""
if not inplace: tods=tods.copy()
todops.fit_basis(tods, dark, highpass=nmode, cuts=cuts, clean_tod=True)
if deslope: utils.deslope(tods, w=8, inplace=True)
return tods
def deproject_vecs_smooth(tods, dark, nmode=50, cuts=None, deslope=True, inplace=True):
if not inplace: tods=tods.copy()
dark = dark.copy()
ftod = fft.rfft(tods)
fdark = fft.rfft(dark)
fdark = todops.smooth_basis_fourier(ftod, fdark)
smooth= np.zeros((fdark.shape[0],dark.shape[1]),dtype=dark.dtype)
fft.ifft(fdark, smooth, normalize=True)
todops.fit_basis(tods, smooth, highpass=nmode, cuts=cuts, clean_tod=True)
if deslope: utils.deslope(tods, w=8, inplace=True)
def filter_common_blockwise(tods, blocks, cuts=None, niter=None,
deslope=True, inplace=True, weight="auto", nmin=5):
# Loop over and filter each block
#np.savetxt("comtod0.txt", tods[0])
if not inplace: tods = tods.copy()
for bi, block in enumerate(blocks):
if len(block) < nmin: continue
btod = np.ascontiguousarray(tods[block])
bcut = None if cuts is None else cuts[block]
common = todops.fit_common(btod, cuts=bcut, niter=niter, clean_tod=True, weight=weight)
tods[block] = btod
#if 0 in block: np.savetxt("comcom.txt", common)
if deslope: utils.deslope(tods, w=8, inplace=True)
#np.savetxt("comtod1.txt", tods[0])
#1/0
return tods
def filter_phase_blockwise(tods, blocks, az, daz=None, cuts=None, niter=None,
deslope=True, inplace=True, weight="auto"):
"""Given a tod[ndet,nsamp], fit for a common azimuth phase signal
per block in blocks[nblock][dets], and subtract it. The binning size
is given in arc minutes."""
# Loop over and filter each block
#np.savetxt("moo0.txt", tods[0])
if not inplace: tods = tods.copy()
for bi, block in enumerate(blocks):
btod = np.ascontiguousarray(tods[block])
bcut = None if cuts is None else cuts[block]
phase = todops.fit_phase_flat(btod, az, daz=daz, cuts=bcut, niter=niter,
clean_tod=True, weight=weight)
tods[block] = btod
if deslope: utils.deslope(tods, w=8, inplace=True)
#np.savetxt("moo.txt", tods[0])
#1/0
return tods