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