Ignore:
Timestamp:
08/03/10 11:41:13 (14 years ago)
Author:
Malte Marquarding
Message:

Tidy up of imports (now imported from asap.). Also fixed some whitespace/tab issues

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/python/simplelinefinder.py

    r1823 r1826  
    11import math
    2 from asap import asaplog
     2from asap.logging import asaplog
    33
    44class StatCalculator:
    … …  
    77       self.s2=0.
    88       self.cnt=0
    9    
     9
    1010   def mean(self):
    1111       if self.cnt<=0:
    1212          raise RuntimeError, "At least one data point has to be defined"
    1313       return self.s/float(self.cnt)
    14    
     14
    1515   def variance(self):
    1616       if self.cnt<=1:
    1717          raise RuntimeError, "At least two data points has to be defined"
    1818       return math.sqrt((self.s2/self.cnt)-(self.s/self.cnt)**2+1e-12)
    19    
     19
    2020   def rms(self):
    2121       """
    … …  
    2525          raise RuntimeError, "At least one data point has to be defined"
    2626       return math.sqrt(self.s2/self.cnt)
    27  
     27
    2828   def add(self, pt):
    2929       self.s = self.s + pt
    … …  
    3636       A simplified class to search for spectral features. The algorithm assumes that the bandpass
    3737       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. 
     38       It works with a list or tuple rather than a scantable and returns the channel pairs.
    3939       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        
     40       it should be used with caution. This class is largely intended to be used with scripts.
     41
    4242       The fully featured version of the algorithm working with scantables is called linefinder.
    4343   '''
    44    
     44
    4545   def __init__(self):
    4646      '''
    47          Initialize the class. 
     47         Initialize the class.
    4848      '''
    4949      self._median = None
    5050      self._rms = None
    51    
     51
    5252   def writeLog(self, str):
    5353      """
    5454         Write user defined string into log file
    55       """ 
     55      """
    5656      asaplog.push(str)
    5757
    5858   def invertChannelSelection(self, nchan, chans, edge = (0,0)):
    5959      """
    60          This method converts a tuple with channel ranges to a tuple which covers all channels 
     60         This method converts a tuple with channel ranges to a tuple which covers all channels
    6161         not selected by the original tuple (optionally edge channels can be discarded)
    62  
    63          nchan - number of channels in the spectrum. 
     62
     63         nchan - number of channels in the spectrum.
    6464         chans - tuple (with even number of elements) containing start and stop channel for all selected ranges
    6565         edge - one number or two element tuple (separate values for two ends) defining how many channels to reject
    … …  
    6969      """
    7070      if nchan<=1:
    71          raise RuntimeError, "Number of channels is supposed to be at least 2, you have %i"% nchan     
     71         raise RuntimeError, "Number of channels is supposed to be at least 2, you have %i"% nchan
    7272      if len(chans)%2!=0:
    7373         raise RuntimeError, "chans is supposed to be a tuple with even number of elements"
    … …  
    106106         spc - tuple with the spectrum (first is the vector of abcissa values, second is the spectrum itself)
    107107         vel_range - a 2-element tuple with start and stop velocity of the range
    108        
     108
    109109         return: a 2-element tuple with channels
    110110         Note, if supplied range is completely outside the spectrum, an empty tuple will be returned
    … …  
    145145      """
    146146         A helper method to compare x and y by absolute value (to do sorting)
    147          
     147
    148148         x - first value
    149149         y - second value
    150          
     150
    151151         return -1,0 or 1 depending on the result of comparison
    152152      """
    … …  
    157157      else:
    158158         return 0
    159          
     159
    160160   def rms(self):
    161161      """
    … …  
    167167         raise RuntimeError, "call find_lines before using the rms method"
    168168      return self._rms
    169    
     169
    170170   def median(self):
    171171      """
    … …  
    176176         raise RuntimeError, "call find_lines before using the median method"
    177177      return self._median
    178    
     178
    179179   def _mergeIntervals(self, lines, spc):
    180180      """
    181181         A helper method to merge intervals.
    182          
     182
    183183         lines - list of tuples with first and last channels of all intervals
    184184         spc - spectrum (to be able to test whether adjacent intervals have the
    … …  
    196196             raise RuntimeError, "this shouldn't have happened!"
    197197          lines.pop(toberemoved[i])
    198          
     198
    199199   def _splitIntervals(self,lines,spc,threshold=3,minchan=3):
    200200      """
    … …  
    203203         Noise is dealt with by taking into account only those extrema, where a difference with
    204204         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.     
     205         class, so the main body of the line detection should be executed first) and there are
     206         at least minchan such points.
    207207      """
    208208      if minchan<1:
    … …  
    219219              if wasIncreasing != None:
    220220                 if isIncreasing != wasIncreasing:
    221                     derivSignReversals.append((ch,isIncreasing)) 
     221                    derivSignReversals.append((ch,isIncreasing))
    222222              wasIncreasing = isIncreasing
    223223          if len(derivSignReversals)==0:
    … …  
    237237             newlines.append((startchan,line[1]))
    238238      return newlines
    239                  
     239
    240240   def find_lines(self,spc,threshold=3,edge=0,minchan=3, tailsearch = True, splitFeatures = False):
    241241      """
    … …  
    243243         taken out perfectly and no channels are flagged within the spectrum. A detection
    244244         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          
     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
    248248         spc - a list or tuple with the spectrum, no default
    249249         threshold - detection threshold, default is 3 sigma, see above for the definition
    … …  
    254254         tailsearch - if True (default), the algorithm attempts to widen each line until
    255255                    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                   
     256                    option off if you need to split the lines according to some criterion
    257257         splitFeatures - if True, the algorithm attempts to split each detected spectral feature into
    258258                    a number of spectral lines (just one local extremum). The default action is
    259259                    not to do it (may give an adverse results if the noise is high)
    260          
     260
    261261         This method returns a list of tuples each containing start and stop 0-based channel
    262262         number of every detected line. Empty list if nothing has been detected.
    263      
     263
    264264         Note. The median and rms about this median is stored inside this class and can
    265265         be obtained with rms and median methods.
    … …  
    270270         raise RuntimeError, "edge is too high (%i), you rejected all channels (%i)" % (edge, len(spc))
    271271      if threshold<=0:
    272          raise RuntimeError, "threshold parameter of find_lines should be positive, you have %s" % threshold     
     272         raise RuntimeError, "threshold parameter of find_lines should be positive, you have %s" % threshold
    273273      if minchan<=0:
    274274         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 
     275
     276      # temporary storage to get statistics, apply edge rejection here
    277277      tmpspc = spc[edge:len(spc)-edge+1]
    278278      if len(tmpspc)<2:
    279          raise RuntimeError, "Too many channels are rejected. Decrease edge parameter or provide a longer spectrum."       
     279         raise RuntimeError, "Too many channels are rejected. Decrease edge parameter or provide a longer spectrum."
    280280      tmpspc.sort()
    281281      self._median=tmpspc[len(tmpspc)/2]
    … …  
    286286      sc = StatCalculator()
    287287      for i in tmpspc[:-int(0.2*len(tmpspc))]:
    288           sc.add(i)         
     288          sc.add(i)
    289289      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      
     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
    294294      #actual line detection
    295295      lines=[]
    … …  
    319319      if nchan>=minchan:
    320320         lines.append((startchan,len(spc)-edge-1))
    321      
     321
    322322      if tailsearch:
    323323         for i in range(len(lines)):
    … …  
    325325             curRange = list(lines[i])
    326326             for x in range(curRange[0],edge,-1):
    327                  isAbove=(spc[x] > self._median)             
     327                 isAbove=(spc[x] > self._median)
    328328                 if wasAbove == None:
    329329                    wasAbove = isAbove
    … …  
    332332                    break
    333333             for x in range(curRange[1],len(spc)-edge):
    334                  isAbove=(spc[x] > self._median)             
     334                 isAbove=(spc[x] > self._median)
    335335                 if isAbove!=wasAbove:
    336336                    curRange[1]=x-1
    337337                    break
    338338             lines[i]=tuple(curRange)
    339          self._mergeIntervals(lines,spc)                           
     339         self._mergeIntervals(lines,spc)
    340340      if splitFeatures:
    341341         return self._splitIntervals(lines,spc,threshold,minchan)
    342       return lines         
     342      return lines
Note: See TracChangeset for help on using the changeset viewer.