source: branches/hpc34/python/asaplinefind.py@ 2708

Last change on this file since 2708 was 2012, checked in by WataruKawasaki, 14 years ago

New Development: Yes

JIRA Issue: Yes (CAS-2373, CAS-2620)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: For Scantable::polyBaseline(), parameters and return value have been changed.

Test Programs:

Put in Release Notes: Yes

Module(s): sdbaseline, sd.linefinder

Description: (1) CAS-2373-related:

(1.1) Modified Scantable::polyBaseline() to have the row-based loop inside.

Now it fits and subtracts baseline for all rows and also output info
about the fitting result to logger and text file, while in the
previous version this method just did baseline fitting/subtraction
for one row only and had to be put inside a row-based loop at the
python side ("poly_baseline()" in scantable.py) and result output had
also to be controlled at the python side. Using a test data from NRO
45m telescope (348,000 rows, 512 channels), the processing time of
scantable.poly_baseline() has reduced from 130 minutes to 5-6 minutes.

(1.2) For accelerating the another polynomial fitting function, namely

scantable.auto_poly_baseline(), added a method
Scantable::autoPolyBaseline(). This basically does the same thing
with Scantable::polyBaseline(), but uses linefinfer also to
automatically flag the line regions.

(1.3) To make linefinder usable in Scantable class, added a method

linefinder.setdata(). This makes it possible to apply linefinder
for a float-list data given as a parameter, without setting scantable,
while it was indispensable to set scantable to use linefinger previously.

(2) CAS-2620-related:

Added Scantable::cubicSplineBaseline() and autoCubicSplineBaseline() for
fit baseline using the cubic spline function. Parameters include npiece
(number of polynomial pieces), clipthresh (clipping threshold), and
clipniter (maximum iteration number).



  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.3 KB
Line 
1import _asap
2
3class linefinder:
4 """
5 The class for automated spectral line search in ASAP.
6
7 Example:
8 fl=linefinder()
9 fl.set_scan(sc)
10 fl.set_options(threshold=3)
11 nlines=fl.find_lines(edge=(50,0))
12 if nlines!=0:
13 print "Found ",nlines," spectral lines"
14 print fl.get_ranges(False)
15 else:
16 print "No lines found!"
17 sc2=sc.poly_baseline(fl.get_mask(),7)
18
19 The algorithm involves a simple threshold criterion. The line is
20 considered to be detected if a specified number of consequtive
21 channels (default is 3) is brighter (with respect to the current baseline
22 estimate) than the threshold times the noise level. This criterion is
23 applied in the iterative procedure updating baseline estimate and trying
24 reduced spectral resolutions to detect broad lines as well. The off-line
25 noise level is determined at each iteration as an average of 80% of the
26 lowest variances across the spectrum (i.e. histogram equalization is
27 used to avoid missing weak lines if strong ones are present). For
28 bad baseline shapes it is recommended to increase the threshold and
29 possibly switch the averaging option off (see set_options) to
30 detect strong lines only, fit a high order baseline and repeat the line
31 search.
32
33 """
34
35 def __init__(self):
36 """
37 Create a line finder object.
38 """
39 self.finder = _asap.linefinder()
40 return
41
42 def set_options(self,threshold=1.7320508075688772,min_nchan=3,
43 avg_limit=8,box_size=0.2,noise_box='all',noise_stat='mean80'):
44 """
45 Set the parameters of the algorithm
46 Parameters:
47 threshold a single channel S/N ratio above which the
48 channel is considered to be a detection
49 Default is sqrt(3), which together with
50 min_nchan=3 gives a 3-sigma criterion
51 min_nchan a minimal number of consequtive channels,
52 which should satisfy a threshold criterion to
53 be a detection. Default is 3.
54 avg_limit A number of consequtive channels not greater than
55 this parameter can be averaged to search for
56 broad lines. Default is 8.
57 box_size A running mean/median box size specified as a fraction
58 of the total spectrum length. Default is 1/5
59 noise_box Area of the spectrum used to estimate noise stats
60 Both string values and numbers are allowed
61 Allowed string values:
62 'all' use all the spectrum (default)
63 'box' noise box is the same as running mean/median
64 box
65 Numeric values are defined as a fraction from the
66 spectrum size. Values should be positive.
67 (noise_box == box_size has the same effect as
68 noise_box = 'box')
69 noise_stat Statistics used to estimate noise, allowed values:
70 'mean80' use the 80% of the lowest deviations
71 in the noise box (default)
72 'median' median of deviations in the noise box
73
74 Note: For bad baselines threshold should be increased,
75 and avg_limit decreased (or even switched off completely by
76 setting this parameter to 1) to avoid detecting baseline
77 undulations instead of real lines.
78 """
79 if noise_stat.lower() not in ["mean80",'median']:
80 raise RuntimeError, "noise_stat argument in linefinder.set_options can only be mean80 or median"
81 nStat = (noise_stat.lower() == "median")
82 nBox = -1.
83 if isinstance(noise_box,str):
84 if noise_box.lower() not in ['all','box']:
85 raise RuntimeError, "string-valued noise_box in linefinder.set_options can only be all or box"
86 if noise_box.lower() == 'box':
87 nBox = box_size
88 else:
89 nBox = float(noise_box)
90 self.finder.setoptions(threshold,min_nchan,avg_limit,box_size,nBox,nStat)
91 return
92
93 def set_scan(self, scan):
94 """
95 Set the 'data' (scantable) to work with.
96 Parameters:
97 scan: a scantable
98 """
99 if not scan:
100 raise RuntimeError, 'Please give a correct scan'
101 self.finder.setscan(scan)
102
103 def set_data(self, spectrum):
104 """
105 Set the 'data' (spectrum) to work with
106 Parameters: a method to allow linefinder work without setting scantable
107 for the purpose of using linefinder inside some method in scantable
108 class. (Dec 22, 2010 by W.Kawasaki)
109 """
110 if isinstance(spectrum, list) or isinstance(spectrum, tuple):
111 if not isinstance(spectrum[0], float):
112 raise RuntimeError, "Parameter 'spectrum' has to be a list or a tuple of float"
113 else:
114 raise RuntimeError, "Parameter 'spectrum' has to be a list or a tuple of float"
115 self.finder.setdata(spectrum)
116
117 def find_lines(self,nRow=0,mask=[],edge=(0,0)):
118 """
119 Search for spectral lines in the scan assigned in set_scan.
120 Parameters:
121 nRow: a row in the scantable to work with
122 mask: an optional mask (e.g. retreived from scantable)
123 edge: an optional number of channels to drop at
124 the edge of the spectrum. If only one value is
125 specified, the same number will be dropped from
126 both sides of the spectrum. Default is to keep
127 all channels
128 A number of lines found will be returned
129 """
130 if isinstance(edge,int):
131 edge=(edge,)
132
133 from asap import _is_sequence_or_number as _is_valid
134
135 if not _is_valid(edge, int):
136 raise RuntimeError, "Parameter 'edge' has to be an integer or \
137 a pair of integers specified as a tuple"
138
139 if len(edge)>2:
140 raise RuntimeError, "The edge parameter should have two \
141 or less elements"
142 return self.finder.findlines(mask,list(edge),nRow)
143 def get_mask(self,invert=False):
144 """
145 Get the mask to mask out all lines that have been found (default)
146
147 Parameters:
148 invert if True, only channels belong to lines will be unmasked
149
150 Note: all channels originally masked by the input mask or
151 dropped out by the edge parameter will still be excluded
152 regardless on the invert option
153 """
154 return self.finder.getmask(invert)
155 def get_ranges(self,defunits=True):
156 """
157 Get ranges (start and end channels or velocities) for all spectral
158 lines found.
159
160 Parameters:
161 defunits if True (default), the range will use the same units
162 as set for the scan (e.g. LSR velocity)
163 if False, the range will be expressed in channels
164 """
165 if (defunits):
166 return self.finder.getlineranges()
167 else:
168 return self.finder.getlinerangesinchannels()
Note: See TracBrowser for help on using the repository browser.