1 | import math |
---|
2 | from asap import asaplog |
---|
3 | |
---|
4 | class StatCalculator: |
---|
5 | def __init__(self): |
---|
6 | self.s=0. |
---|
7 | self.s2=0. |
---|
8 | self.cnt=0 |
---|
9 | |
---|
10 | def mean(self): |
---|
11 | if self.cnt<=0: |
---|
12 | raise RuntimeError, "At least one data point has to be defined" |
---|
13 | return self.s/float(self.cnt) |
---|
14 | |
---|
15 | def variance(self): |
---|
16 | if self.cnt<=1: |
---|
17 | raise RuntimeError, "At least two data points has to be defined" |
---|
18 | return math.sqrt((self.s2/self.cnt)-(self.s/self.cnt)**2+1e-12) |
---|
19 | |
---|
20 | def rms(self): |
---|
21 | """ |
---|
22 | return rms of the accumulated sample |
---|
23 | """ |
---|
24 | if self.cnt<=0: |
---|
25 | raise RuntimeError, "At least one data point has to be defined" |
---|
26 | return math.sqrt(self.s2/self.cnt) |
---|
27 | |
---|
28 | def add(self, pt): |
---|
29 | self.s = self.s + pt |
---|
30 | self.s2 = self.s2 + pt*pt |
---|
31 | self.cnt = self.cnt + 1 |
---|
32 | |
---|
33 | |
---|
34 | class simplelinefinder: |
---|
35 | ''' |
---|
36 | A simplified class to search for spectral features. The algorithm assumes that the bandpass |
---|
37 | is taken out perfectly and no spectral channels are flagged (except some edge channels). |
---|
38 | It works with a list or tuple rather than a scantable and returns the channel pairs. |
---|
39 | There is an optional feature to attempt to split the detected lines into components, although |
---|
40 | it should be used with caution. This class is largely intended to be used with scripts. |
---|
41 | |
---|
42 | The fully featured version of the algorithm working with scantables is called linefinder. |
---|
43 | ''' |
---|
44 | |
---|
45 | def __init__(self): |
---|
46 | ''' |
---|
47 | Initialize the class. |
---|
48 | ''' |
---|
49 | self._median = None |
---|
50 | self._rms = None |
---|
51 | |
---|
52 | def writeLog(self, str): |
---|
53 | """ |
---|
54 | Write user defined string into log file |
---|
55 | """ |
---|
56 | asaplog.push(str) |
---|
57 | |
---|
58 | def invertChannelSelection(self, nchan, chans, edge = (0,0)): |
---|
59 | """ |
---|
60 | This method converts a tuple with channel ranges to a tuple which covers all channels |
---|
61 | not selected by the original tuple (optionally edge channels can be discarded) |
---|
62 | |
---|
63 | nchan - number of channels in the spectrum. |
---|
64 | chans - tuple (with even number of elements) containing start and stop channel for all selected ranges |
---|
65 | edge - one number or two element tuple (separate values for two ends) defining how many channels to reject |
---|
66 | |
---|
67 | return: a tuple with inverted channel selection |
---|
68 | Note, at this stage channel ranges are assumed to be sorted and without overlap |
---|
69 | """ |
---|
70 | if nchan<=1: |
---|
71 | raise RuntimeError, "Number of channels is supposed to be at least 2, you have %i"% nchan |
---|
72 | if len(chans)%2!=0: |
---|
73 | raise RuntimeError, "chans is supposed to be a tuple with even number of elements" |
---|
74 | tempedge = edge |
---|
75 | if not isinstance(tempedge,tuple): |
---|
76 | tempedge = (edge,edge) |
---|
77 | if len(tempedge)!=2: |
---|
78 | raise RuntimeError, "edge parameter is supposed to be a two-element tuple or a single number" |
---|
79 | if tempedge[0]<0 or tempedge[1]<0: |
---|
80 | raise RuntimeError, "number of edge rejected cannels is supposed to be positive" |
---|
81 | startchan = tempedge[0] |
---|
82 | stopchan = nchan - tempedge[1] |
---|
83 | if stopchan-startchan<0: |
---|
84 | raise RuntimeError, "Edge rejection rejected all channels" |
---|
85 | ranges = [] |
---|
86 | curstart = startchan |
---|
87 | for i in range(0,len(chans),2): |
---|
88 | if chans[i+1]<curstart: |
---|
89 | continue |
---|
90 | elif chans[i]<=curstart: |
---|
91 | curstart = chans[i+1]+1 |
---|
92 | else: |
---|
93 | ranges.append(curstart) |
---|
94 | ranges.append(chans[i]-1) |
---|
95 | curstart = chans[i+1]+1 |
---|
96 | if curstart<stopchan: |
---|
97 | ranges.append(curstart) |
---|
98 | ranges.append(stopchan-1) |
---|
99 | return tuple(ranges) |
---|
100 | |
---|
101 | def channelRange(self, spc, vel_range): |
---|
102 | """ |
---|
103 | A helper method which works on a tuple with abcissa/flux vectors (i.e. as returned by uvSpectrum). It |
---|
104 | allows to convert supplied velocity range into the channel range. |
---|
105 | |
---|
106 | spc - tuple with the spectrum (first is the vector of abcissa values, second is the spectrum itself) |
---|
107 | vel_range - a 2-element tuple with start and stop velocity of the range |
---|
108 | |
---|
109 | return: a 2-element tuple with channels |
---|
110 | Note, if supplied range is completely outside the spectrum, an empty tuple will be returned |
---|
111 | """ |
---|
112 | if len(spc) != 2: |
---|
113 | raise RuntimeError, "spc supplied to channelRange is supposed to have 2 elements" |
---|
114 | if len(spc[0]) != len(spc[1]): |
---|
115 | raise RuntimeError, "spc supplied to channelRange is supposed to have 2 vectors of equal length" |
---|
116 | if len(vel_range) != 2: |
---|
117 | raise RuntimeError, "vel_range supplied to channelRange is supposed to have 2 elements" |
---|
118 | if vel_range[0]>=vel_range[1]: |
---|
119 | raise RuntimeError, "start velocity is supposed to be less than end velocity, vel_range: %s" % vel_range |
---|
120 | if len(spc[0])<=2: |
---|
121 | raise RuntimeError, "Spectrum should contain more than 2 points, you have %i" % len(spc[0]) |
---|
122 | chans = list(vel_range) |
---|
123 | for j in range(len(chans)): |
---|
124 | chans[j] = -1 |
---|
125 | for i in range(len(spc[0])): |
---|
126 | if i!=0: |
---|
127 | prev_vel = spc[0][i-1] |
---|
128 | else: |
---|
129 | prev_vel = spc[0][i+1] |
---|
130 | delta = max(prev_vel, spc[0][i]) - min(prev_vel, spc[0][i]) |
---|
131 | for j in range(len(vel_range)): |
---|
132 | if abs(vel_range[j]-spc[0][i])<delta: |
---|
133 | chans[j] = i |
---|
134 | if chans[0] == chans[1]: |
---|
135 | return () |
---|
136 | if chans[1] == -1: |
---|
137 | chans[1] = len(spc[0])-1 |
---|
138 | if chans[0] == -1: |
---|
139 | chans[0] = 0 |
---|
140 | if chans[0]<=chans[1]: |
---|
141 | return tuple(chans) |
---|
142 | return (chans[1],chans[0]) |
---|
143 | |
---|
144 | def _absvalComp(self,x,y): |
---|
145 | """ |
---|
146 | A helper method to compare x and y by absolute value (to do sorting) |
---|
147 | |
---|
148 | x - first value |
---|
149 | y - second value |
---|
150 | |
---|
151 | return -1,0 or 1 depending on the result of comparison |
---|
152 | """ |
---|
153 | if abs(x)<abs(y): |
---|
154 | return -1 |
---|
155 | elif abs(x)>abs(y): |
---|
156 | return 1 |
---|
157 | else: |
---|
158 | return 0 |
---|
159 | |
---|
160 | def rms(self): |
---|
161 | """ |
---|
162 | Return rms scatter of the spectral points (with respect to the median) calculated |
---|
163 | during last find_lines call. Note, this method throws an exception if |
---|
164 | find_lines has never been called. |
---|
165 | """ |
---|
166 | if self._rms==None: |
---|
167 | raise RuntimeError, "call find_lines before using the rms method" |
---|
168 | return self._rms |
---|
169 | |
---|
170 | def median(self): |
---|
171 | """ |
---|
172 | Return the median of the last spectrum passed to find_lines. |
---|
173 | Note, this method throws an exception if find_lines has never been called. |
---|
174 | """ |
---|
175 | if self._median==None: |
---|
176 | raise RuntimeError, "call find_lines before using the median method" |
---|
177 | return self._median |
---|
178 | |
---|
179 | def _mergeIntervals(self, lines, spc): |
---|
180 | """ |
---|
181 | A helper method to merge intervals. |
---|
182 | |
---|
183 | lines - list of tuples with first and last channels of all intervals |
---|
184 | spc - spectrum (to be able to test whether adjacent intervals have the |
---|
185 | same sign. |
---|
186 | """ |
---|
187 | toberemoved = [] |
---|
188 | for i in range(1,len(lines)): |
---|
189 | if lines[i-1][1]+1>=lines[i][0]: |
---|
190 | if (spc[lines[i-1][1]]>self._median) == (spc[lines[i][0]]>self._median): |
---|
191 | lines[i] = (lines[i-1][0],lines[i][1]) |
---|
192 | toberemoved.append(i-1) |
---|
193 | toberemoved.sort() |
---|
194 | for i in range(len(toberemoved)-1,-1,-1): |
---|
195 | if toberemoved[i]>=len(lines): |
---|
196 | raise RuntimeError, "this shouldn't have happened!" |
---|
197 | lines.pop(toberemoved[i]) |
---|
198 | |
---|
199 | def _splitIntervals(self,lines,spc,threshold=3,minchan=3): |
---|
200 | """ |
---|
201 | A helper method used in the spectral line detection routine. It splits given intervals |
---|
202 | into a number of "spectral lines". Each line is characterised by a single extremum. |
---|
203 | Noise is dealt with by taking into account only those extrema, where a difference with |
---|
204 | respect to surrounding spectral points exceeds threshold times rms (stored inside this |
---|
205 | class, so the main body of the line detection should be executed first) and there are |
---|
206 | at least minchan such points. |
---|
207 | """ |
---|
208 | if minchan<1: |
---|
209 | raise RuntimeError, "minchan in _splitIntervals is not supposed to be less than 1, you have %s" % minchan |
---|
210 | newlines = [] |
---|
211 | for line in lines: |
---|
212 | if line[1]-line[0]+1 <= minchan: |
---|
213 | newlines.append(line) |
---|
214 | wasIncreasing = None |
---|
215 | derivSignReversals = [] |
---|
216 | for ch in range(line[0]+1,line[1]+1): |
---|
217 | diff=spc[ch]-spc[ch-1] |
---|
218 | isIncreasing = (diff>0) |
---|
219 | if wasIncreasing != None: |
---|
220 | if isIncreasing != wasIncreasing: |
---|
221 | derivSignReversals.append((ch,isIncreasing)) |
---|
222 | wasIncreasing = isIncreasing |
---|
223 | if len(derivSignReversals)==0: |
---|
224 | newlines.append(line) |
---|
225 | elif len(derivSignReversals)%2 != 1: |
---|
226 | self.writeLog("SplitIntervals warning. Derivative is expected to have odd number of reversals within the interval: \"%s\" " % derivSignReversals); |
---|
227 | newlines.append(line) |
---|
228 | elif derivSignReversals[0][1]!=derivSignReversals[-1][1]: |
---|
229 | self.writeLog("SplitIntervals warning. Derivative is expected to have the same sign at the start and at the end of each interval: \"%s\"" % derivSignReversals) |
---|
230 | newlines.append(line) |
---|
231 | else: |
---|
232 | startchan = line[0] |
---|
233 | for i in range(len(derivSignReversals)): |
---|
234 | if i%2 == 1: |
---|
235 | newlines.append((startchan,derivSignReversals[i][0]-1)) |
---|
236 | startchan = derivSignReversals[i][0] |
---|
237 | newlines.append((startchan,line[1])) |
---|
238 | return newlines |
---|
239 | |
---|
240 | def find_lines(self,spc,threshold=3,edge=0,minchan=3, tailsearch = True, splitFeatures = False): |
---|
241 | """ |
---|
242 | A simple spectral line detection routine, which assumes that bandpass has been |
---|
243 | taken out perfectly and no channels are flagged within the spectrum. A detection |
---|
244 | is reported if consequtive minchan number of channels is consistently above or |
---|
245 | below the median value. The threshold is given in terms of the rms calculated |
---|
246 | using 80% of the lowest data points by the absolute value (with respect to median) |
---|
247 | |
---|
248 | spc - a list or tuple with the spectrum, no default |
---|
249 | threshold - detection threshold, default is 3 sigma, see above for the definition |
---|
250 | edge - if non-zero, this number of spectral channels will be rejected at the edge. |
---|
251 | Default is not to do any rejection. |
---|
252 | minchan - minimum number of consequitive channels exceeding threshold to claim the |
---|
253 | detection, default is 3. |
---|
254 | tailsearch - if True (default), the algorithm attempts to widen each line until |
---|
255 | its flux crosses the median. It merges lines if necessary. Set this |
---|
256 | option off if you need to split the lines according to some criterion |
---|
257 | splitFeatures - if True, the algorithm attempts to split each detected spectral feature into |
---|
258 | a number of spectral lines (just one local extremum). The default action is |
---|
259 | not to do it (may give an adverse results if the noise is high) |
---|
260 | |
---|
261 | This method returns a list of tuples each containing start and stop 0-based channel |
---|
262 | number of every detected line. Empty list if nothing has been detected. |
---|
263 | |
---|
264 | Note. The median and rms about this median is stored inside this class and can |
---|
265 | be obtained with rms and median methods. |
---|
266 | """ |
---|
267 | if edge<0: |
---|
268 | raise RuntimeError, "edge parameter of find_lines should be non-negative, you have %s" % edge |
---|
269 | if 2*edge>=len(spc): |
---|
270 | raise RuntimeError, "edge is too high (%i), you rejected all channels (%i)" % (edge, len(spc)) |
---|
271 | if threshold<=0: |
---|
272 | raise RuntimeError, "threshold parameter of find_lines should be positive, you have %s" % threshold |
---|
273 | if minchan<=0: |
---|
274 | raise RuntimeError, "minchan parameter of find_lines should be positive, you have %s" % minchan |
---|
275 | |
---|
276 | # temporary storage to get statistics, apply edge rejection here |
---|
277 | tmpspc = spc[edge:len(spc)-edge+1] |
---|
278 | if len(tmpspc)<2: |
---|
279 | raise RuntimeError, "Too many channels are rejected. Decrease edge parameter or provide a longer spectrum." |
---|
280 | tmpspc.sort() |
---|
281 | self._median=tmpspc[len(tmpspc)/2] |
---|
282 | # work with offsets from the median and sort by absolute values |
---|
283 | for i in range(len(tmpspc)): |
---|
284 | tmpspc[i]-=self._median |
---|
285 | tmpspc.sort(cmp=self._absvalComp) |
---|
286 | sc = StatCalculator() |
---|
287 | for i in tmpspc[:-int(0.2*len(tmpspc))]: |
---|
288 | sc.add(i) |
---|
289 | self._rms=sc.rms() |
---|
290 | |
---|
291 | self.writeLog("Spectral line detection with edge=%i, threshold=%f, minchan=%i and tailsearch=%s" % (edge,threshold, minchan, tailsearch)) |
---|
292 | self.writeLog("statistics: median=%f, rms=%f" % (self._median, self._rms)) |
---|
293 | |
---|
294 | #actual line detection |
---|
295 | lines=[] |
---|
296 | wasAbove = None |
---|
297 | nchan = 0 |
---|
298 | startchan=None |
---|
299 | for i in range(edge,len(spc)-edge): |
---|
300 | if abs(spc[i]-self._median)>threshold*self._rms: |
---|
301 | isAbove=(spc[i] > self._median) |
---|
302 | if nchan!=0: |
---|
303 | if wasAbove == isAbove: |
---|
304 | nchan+=1 |
---|
305 | else: |
---|
306 | if nchan>=minchan: |
---|
307 | lines.append((startchan,i-1)) |
---|
308 | nchan=1 |
---|
309 | wasAbove = isAbove |
---|
310 | startchan = i |
---|
311 | else: |
---|
312 | nchan=1 |
---|
313 | wasAbove = isAbove |
---|
314 | startchan = i |
---|
315 | else: |
---|
316 | if nchan>=minchan: |
---|
317 | lines.append((startchan,i-1)) |
---|
318 | nchan = 0 |
---|
319 | if nchan>=minchan: |
---|
320 | lines.append((startchan,len(spc)-edge-1)) |
---|
321 | |
---|
322 | if tailsearch: |
---|
323 | for i in range(len(lines)): |
---|
324 | wasAbove = None |
---|
325 | curRange = list(lines[i]) |
---|
326 | for x in range(curRange[0],edge,-1): |
---|
327 | isAbove=(spc[x] > self._median) |
---|
328 | if wasAbove == None: |
---|
329 | wasAbove = isAbove |
---|
330 | if isAbove!=wasAbove: |
---|
331 | curRange[0]=x+1 |
---|
332 | break |
---|
333 | for x in range(curRange[1],len(spc)-edge): |
---|
334 | isAbove=(spc[x] > self._median) |
---|
335 | if isAbove!=wasAbove: |
---|
336 | curRange[1]=x-1 |
---|
337 | break |
---|
338 | lines[i]=tuple(curRange) |
---|
339 | self._mergeIntervals(lines,spc) |
---|
340 | if splitFeatures: |
---|
341 | return self._splitIntervals(lines,spc,threshold,minchan) |
---|
342 | return lines |
---|