Changeset 2186 for trunk/python
- Timestamp:
- 06/07/11 23:49:53 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/python/scantable.py
r2178 r2186 76 76 else: 77 77 return False 78 79 def normalise_edge_param(edge): 80 """\ 81 Convert a given edge value to a one-dimensional array that can be 82 given to baseline-fitting/subtraction functions. 83 The length of the output value will be an even because values for 84 the both sides of spectra are to be contained for each IF. When 85 the length is 2, the values will be applied to all IFs. If the length 86 is larger than 2, it will be 2*ifnos(). 87 Accepted format of edge include: 88 * an integer - will be used for both sides of spectra of all IFs. 89 e.g. 10 is converted to [10,10] 90 * an empty list/tuple [] - converted to [0, 0] and used for all IFs. 91 * a list/tuple containing an integer - same as the above case. 92 e.g. [10] is converted to [10,10] 93 * a list/tuple containing two integers - will be used for all IFs. 94 e.g. [5,10] is output as it is. no need to convert. 95 * a list/tuple of lists/tuples containing TWO integers - 96 each element of edge will be used for each IF. 97 e.g. [[5,10],[15,20]] - [5,10] for IF[0] and [15,20] for IF[1]. 98 99 If an element contains the same integer values, the input 'edge' 100 parameter can be given in a simpler shape in the following cases: 101 ** when len(edge)!=2 102 any elements containing the same values can be replaced 103 to single integers. 104 e.g. [[15,15]] can be simplified to [15] (or [15,15] or 15 also). 105 e.g. [[1,1],[2,2],[3,3]] can be simplified to [1,2,3]. 106 ** when len(edge)=2 107 care is needed for this case: ONLY ONE of the 108 elements can be a single integer, 109 e.g. [[5,5],[10,10]] can be simplified to [5,[10,10]] 110 or [[5,5],10], but can NOT be simplified to [5,10]. 111 when [5,10] given, it is interpreted as 112 [[5,10],[5,10],[5,10],...] instead, as shown before. 113 """ 114 from asap import _is_sequence_or_number as _is_valid 115 if isinstance(edge, list) or isinstance(edge, tuple): 116 for edgepar in edge: 117 if not _is_valid(edgepar, int): 118 raise ValueError, "Each element of the 'edge' tuple has \ 119 to be a pair of integers or an integer." 120 if isinstance(edgepar, list) or isinstance(edgepar, tuple): 121 if len(edgepar) != 2: 122 raise ValueError, "Each element of the 'edge' tuple has \ 123 to be a pair of integers or an integer." 124 else: 125 if not _is_valid(edge, int): 126 raise ValueError, "Parameter 'edge' has to be an integer or a \ 127 pair of integers specified as a tuple. \ 128 Nested tuples are allowed \ 129 to make individual selection for different IFs." 130 131 132 if isinstance(edge, int): 133 edge = [ edge, edge ] # e.g. 3 => [3,3] 134 elif isinstance(edge, list) or isinstance(edge, tuple): 135 if len(edge) == 0: 136 edge = [0, 0] # e.g. [] => [0,0] 137 elif len(edge) == 1: 138 if isinstance(edge[0], int): 139 edge = [ edge[0], edge[0] ] # e.g. [1] => [1,1] 140 141 commonedge = True 142 if len(edge) > 2: commonedge = False 143 else: 144 for edgepar in edge: 145 if isinstance(edgepar, list) or isinstance(edgepar, tuple): 146 commonedge = False 147 break 148 149 if commonedge: 150 if len(edge) > 1: 151 norm_edge = edge 152 else: 153 norm_edge = edge + edge 154 else: 155 norm_edge = [] 156 for edgepar in edge: 157 if isinstance(edgepar, int): 158 norm_edge += [edgepar, edgepar] 159 else: 160 norm_edge += edgepar 161 162 return norm_edge 163 164 def raise_fitting_failure_exception(e): 165 msg = "The fit failed, possibly because it didn't converge." 166 if rcParams["verbose"]: 167 asaplog.push(str(e)) 168 asaplog.push(str(msg)) 169 else: 170 raise RuntimeError(str(e)+'\n'+msg) 78 171 172 79 173 class scantable(Scantable): 80 174 """\ … … 1120 1214 1121 1215 @asaplog_post_dec 1122 def fft(self, getrealimag=False, rowno=[]):1216 def fft(self, rowno=[], mask=[], getrealimag=False): 1123 1217 """\ 1124 1218 Apply FFT to the spectra. … … 1126 1220 1127 1221 Parameters: 1222 1223 rowno: The row number(s) to be processed. int, list 1224 and tuple are accepted. By default ([]), FFT 1225 is applied to the whole data. 1226 1227 mask: Auxiliary channel mask(s). Given as a boolean 1228 list, it is applied to all specified rows. 1229 A list of boolean lists can also be used to 1230 apply different masks. In the latter case, the 1231 length of 'mask' must be the same as that of 1232 'rowno'. The default is []. 1128 1233 1129 1234 getrealimag: If True, returns the real and imaginary part … … 1133 1238 phase (argument, in unit of radian). 1134 1239 1135 rowno: The row number(s) to be processed. int, list1136 and tuple are accepted. By default ([]), FFT1137 is applied to the whole data.1138 1139 1240 Returns: 1140 1241 1141 A dictionary containing two keys, i.e., 'real' and 'imag' for 1142 getrealimag = True, or 'ampl' and 'phase' for getrealimag = False, 1242 A list of dictionaries containing the results for each spectrum. 1243 Each dictionary contains two values, the real and the imaginary 1244 parts when getrealimag = True, or the amplitude(absolute value) 1245 and the phase(argument) when getrealimag = False. The key for 1246 these values are 'real' and 'imag', or 'ampl' and 'phase', 1143 1247 respectively. 1144 The value for each key is a list of lists containing the FFT 1145 results from each spectrum. 1146 1147 """ 1148 if getrealimag: 1149 res = {"real":[], "imag":[]} # return real and imaginary values 1150 else: 1151 res = {"ampl":[], "phase":[]} # return amplitude and phase(argument) 1152 1248 """ 1153 1249 if isinstance(rowno, int): 1154 1250 rowno = [rowno] 1155 1251 elif not (isinstance(rowno, list) or isinstance(rowno, tuple)): 1156 raise TypeError("The row number(s) should be int, list or tuple.") 1252 raise TypeError("The row number(s) must be int, list or tuple.") 1253 1254 if len(rowno) == 0: rowno = [i for i in xrange(self.nrow())] 1255 1256 if not (isinstance(mask, list) or isinstance(mask, tuple)): 1257 raise TypeError("The mask must be a boolean list or a list of boolean list.") 1258 if len(mask) == 0: mask = [True for i in xrange(self.nchan())] 1259 if isinstance(mask[0], bool): mask = [mask] 1260 elif not (isinstance(mask[0], list) or isinstance(mask[0], tuple)): 1261 raise TypeError("The mask must be a boolean list or a list of boolean list.") 1262 1263 usecommonmask = (len(mask) == 1) 1264 if not usecommonmask: 1265 if len(mask) != len(rowno): 1266 raise ValueError("When specifying masks for each spectrum, the numbers of them must be identical.") 1267 for amask in mask: 1268 if len(amask) != self.nchan(): 1269 raise ValueError("The spectra and the mask have different length.") 1157 1270 1158 nrow = len(rowno) 1159 if nrow == 0: nrow = self.nrow() # by default, calculate for all rows. 1160 1161 fspec = self._math._fft(self, rowno, getrealimag) 1162 nspec = len(fspec)/nrow 1163 1164 i = 0 1165 while (i < nrow): 1166 v1 = [] 1167 v2 = [] 1271 res = [] 1272 1273 imask = 0 1274 for whichrow in rowno: 1275 fspec = self._fft(whichrow, mask[imask], getrealimag) 1276 nspec = len(fspec) 1168 1277 1169 j = 0 1170 while (j < nspec): 1171 k = i*nspec + j 1172 v1.append(fspec[k]) 1173 v2.append(fspec[k+1]) 1174 j += 2 1278 i=0 1279 v1=[] 1280 v2=[] 1281 reselem = {"real":[],"imag":[]} if getrealimag else {"ampl":[],"phase":[]} 1282 1283 while (i < nspec): 1284 v1.append(fspec[i]) 1285 v2.append(fspec[i+1]) 1286 i+=2 1175 1287 1176 1288 if getrealimag: 1177 res ["real"].append(v1)1178 res ["imag"].append(v2)1289 reselem["real"] += v1 1290 reselem["imag"] += v2 1179 1291 else: 1180 res ["ampl"].append(v1)1181 res ["phase"].append(v2)1292 reselem["ampl"] += v1 1293 reselem["phase"] += v2 1182 1294 1183 i += 1 1184 1295 res.append(reselem) 1296 1297 if not usecommonmask: imask += 1 1298 1185 1299 return res 1186 1300 … … 2131 2245 else: return s 2132 2246 2133 2134 @asaplog_post_dec 2135 def sinusoid_baseline(self, insitu=None, mask=None, nwave=None, maxwavelength=None, 2136 clipthresh=None, clipniter=None, plot=None, getresidual=None, outlog=None, blfile=None): 2247 @asaplog_post_dec 2248 def _parse_wn(self, wn): 2249 if isinstance(wn, list) or isinstance(wn, tuple): 2250 return wn 2251 elif isinstance(wn, int): 2252 return [ wn ] 2253 elif isinstance(wn, str): 2254 if '-' in wn: # case 'a-b' : return [a,a+1,...,b-1,b] 2255 val = wn.split('-') 2256 val = [int(val[0]), int(val[1])] 2257 val.sort() 2258 res = [i for i in xrange(val[0], val[1]+1)] 2259 elif wn[:2] == '<=' or wn[:2] == '=<': # cases '<=a','=<a' : return [0,1,...,a-1,a] 2260 val = int(wn[2:])+1 2261 res = [i for i in xrange(val)] 2262 elif wn[-2:] == '>=' or wn[-2:] == '=>': # cases 'a>=','a=>' : return [0,1,...,a-1,a] 2263 val = int(wn[:-2])+1 2264 res = [i for i in xrange(val)] 2265 elif wn[0] == '<': # case '<a' : return [0,1,...,a-2,a-1] 2266 val = int(wn[1:]) 2267 res = [i for i in xrange(val)] 2268 elif wn[-1] == '>': # case 'a>' : return [0,1,...,a-2,a-1] 2269 val = int(wn[:-1]) 2270 res = [i for i in xrange(val)] 2271 elif wn[:2] == '>=' or wn[:2] == '=>': # cases '>=a','=>a' : return [a,a+1,...,a_nyq] 2272 val = int(wn[2:]) 2273 res = [i for i in xrange(val, self.nchan()/2+1)] 2274 elif wn[-2:] == '<=' or wn[-2:] == '=<': # cases 'a<=','a=<' : return [a,a+1,...,a_nyq] 2275 val = int(wn[:-2]) 2276 res = [i for i in xrange(val, self.nchan()/2+1)] 2277 elif wn[0] == '>': # case '>a' : return [a+1,a+2,...,a_nyq] 2278 val = int(wn[1:])+1 2279 res = [i for i in xrange(val, self.nchan()/2+1)] 2280 elif wn[-1] == '<': # case 'a<' : return [a+1,a+2,...,a_nyq] 2281 val = int(wn[:-1])+1 2282 res = [i for i in xrange(val, self.nchan()/2+1)] 2283 2284 return res 2285 else: 2286 msg = 'wrong value given for addwn/rejwn' 2287 raise RuntimeError(msg) 2288 2289 2290 @asaplog_post_dec 2291 def sinusoid_baseline(self, insitu=None, mask=None, applyfft=None, fftmethod=None, fftthresh=None, 2292 addwn=None, rejwn=None, clipthresh=None, clipniter=None, plot=None, 2293 getresidual=None, outlog=None, blfile=None): 2137 2294 """\ 2138 2295 Return a scan which has been baselined (all rows) with sinusoidal functions. 2139 2296 Parameters: 2140 insitu: If False a new scantable is returned.2297 insitu: if False a new scantable is returned. 2141 2298 Otherwise, the scaling is done in-situ 2142 2299 The default is taken from .asaprc (False) 2143 mask: An optional mask 2144 nwave: the maximum wave number of sinusoids within 2145 maxwavelength * (spectral range). 2146 The default is 3 (i.e., sinusoids with wave 2147 number of 0(=constant), 1, 2, and 3 are 2148 used for fitting). Also it is possible to 2149 explicitly specify all the wave numbers to 2150 be used, by giving a list including them 2151 (e.g. [0,1,2,15,16]). 2152 maxwavelength: the longest sinusoidal wavelength. The 2153 default is 1.0 (unit: spectral range). 2300 mask: an optional mask 2301 applyfft: if True use some method, such as FFT, to find 2302 strongest sinusoidal components in the wavenumber 2303 domain to be used for baseline fitting. 2304 default is True. 2305 fftmethod: method to find the strong sinusoidal components. 2306 now only 'fft' is available and it is the default. 2307 fftthresh: the threshold to select wave numbers to be used for 2308 fitting from the distribution of amplitudes in the 2309 wavenumber domain. 2310 both float and string values accepted. 2311 given a float value, the unit is set to sigma. 2312 for string values, allowed formats include: 2313 'xsigma' or 'x' (= x-sigma level. e.g., '3sigma'), or 2314 'topx' (= the x strongest ones, e.g. 'top5'). 2315 default is 3.0 (unit: sigma). 2316 addwn: the additional wave numbers to be used for fitting. 2317 list or integer value is accepted to specify every 2318 wave numbers. also string value can be used in case 2319 you need to specify wave numbers in a certain range, 2320 e.g., 'a-b' (= a, a+1, a+2, ..., b-1, b), 2321 '<a' (= 0,1,...,a-2,a-1), 2322 '>=a' (= a, a+1, ... up to the maximum wave 2323 number corresponding to the Nyquist 2324 frequency for the case of FFT). 2325 default is []. 2326 rejwn: the wave numbers NOT to be used for fitting. 2327 can be set just as addwn but has higher priority: 2328 wave numbers which are specified both in addwn 2329 and rejwn will NOT be used. default is []. 2154 2330 clipthresh: Clipping threshold. (default is 3.0, unit: sigma) 2155 2331 clipniter: maximum number of iteration of 'clipthresh'-sigma clipping (default is 0) … … 2164 2340 blfile: Name of a text file in which the best-fit 2165 2341 parameter values to be written 2166 (default is "": no file/logger output)2342 (default is '': no file/logger output) 2167 2343 2168 2344 Example: … … 2170 2346 # wave numbers in spectral window up to 10, 2171 2347 # also with 3-sigma clipping, iteration up to 4 times 2172 bscan = scan.sinusoid_baseline( nwave=10,clipthresh=3.0,clipniter=4)2348 bscan = scan.sinusoid_baseline(addwn='<=10',clipthresh=3.0,clipniter=4) 2173 2349 2174 2350 Note: … … 2177 2353 """ 2178 2354 2179 varlist = vars() 2355 try: 2356 varlist = vars() 2180 2357 2181 if insitu is None: insitu = rcParams["insitu"] 2182 if insitu: 2183 workscan = self 2184 else: 2185 workscan = self.copy() 2186 2187 nchan = workscan.nchan() 2188 2189 if mask is None: mask = [True for i in xrange(nchan)] 2190 if nwave is None: nwave = 3 2191 if maxwavelength is None: maxwavelength = 1.0 2192 if clipthresh is None: clipthresh = 3.0 2193 if clipniter is None: clipniter = 0 2194 if plot is None: plot = False 2195 if getresidual is None: getresidual = True 2196 if outlog is None: outlog = False 2197 if blfile is None: blfile = "" 2198 2199 if isinstance(nwave, int): 2200 in_nwave = nwave 2201 nwave = [] 2202 for i in xrange(in_nwave+1): nwave.append(i) 2203 2204 outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile))) 2205 2206 try: 2358 if insitu is None: insitu = rcParams['insitu'] 2359 if insitu: 2360 workscan = self 2361 else: 2362 workscan = self.copy() 2363 2364 if mask is None: mask = [True for i in xrange(workscan.nchan())] 2365 if applyfft is None: applyfft = True 2366 if fftmethod is None: fftmethod = 'fft' 2367 if fftthresh is None: fftthresh = 3.0 2368 if addwn is None: addwn = [] 2369 if rejwn is None: rejwn = [] 2370 if clipthresh is None: clipthresh = 3.0 2371 if clipniter is None: clipniter = 0 2372 if plot is None: plot = False 2373 if getresidual is None: getresidual = True 2374 if outlog is None: outlog = False 2375 if blfile is None: blfile = '' 2376 2207 2377 #CURRENTLY, PLOT=true is UNAVAILABLE UNTIL sinusoidal fitting is implemented as a fitter method. 2208 workscan._sinusoid_baseline(mask, nwave, maxwavelength, clipthresh, clipniter, getresidual, outlog, blfile) 2209 2210 workscan._add_history("sinusoid_baseline", varlist) 2378 workscan._sinusoid_baseline(mask, applyfft, fftmethod.lower(), str(fftthresh).lower(), workscan._parse_wn(addwn), workscan._parse_wn(rejwn), clipthresh, clipniter, getresidual, outlog, blfile) 2379 workscan._add_history('sinusoid_baseline', varlist) 2211 2380 2212 2381 if insitu: … … 2216 2385 2217 2386 except RuntimeError, e: 2218 msg = "The fit failed, possibly because it didn't converge." 2219 if rcParams["verbose"]: 2220 asaplog.push(str(e)) 2221 asaplog.push(str(msg)) 2222 return 2223 else: 2224 raise RuntimeError(str(e)+'\n'+msg) 2225 2226 2227 def auto_sinusoid_baseline(self, insitu=None, mask=None, nwave=None, maxwavelength=None, 2228 clipthresh=None, clipniter=None, edge=None, threshold=None, 2387 raise_fitting_failure_exception(e) 2388 2389 2390 @asaplog_post_dec 2391 def auto_sinusoid_baseline(self, insitu=None, mask=None, applyfft=None, fftmethod=None, fftthresh=None, 2392 addwn=None, rejwn=None, clipthresh=None, clipniter=None, edge=None, threshold=None, 2229 2393 chan_avg_limit=None, plot=None, getresidual=None, outlog=None, blfile=None): 2230 2394 """\ … … 2238 2402 The default is taken from .asaprc (False) 2239 2403 mask: an optional mask retreived from scantable 2240 nwave: the maximum wave number of sinusoids within 2241 maxwavelength * (spectral range). 2242 The default is 3 (i.e., sinusoids with wave 2243 number of 0(=constant), 1, 2, and 3 are 2244 used for fitting). Also it is possible to 2245 explicitly specify all the wave numbers to 2246 be used, by giving a list including them 2247 (e.g. [0,1,2,15,16]). 2248 maxwavelength: the longest sinusoidal wavelength. The 2249 default is 1.0 (unit: spectral range). 2404 applyfft: if True use some method, such as FFT, to find 2405 strongest sinusoidal components in the wavenumber 2406 domain to be used for baseline fitting. 2407 default is True. 2408 fftmethod: method to find the strong sinusoidal components. 2409 now only 'fft' is available and it is the default. 2410 fftthresh: the threshold to select wave numbers to be used for 2411 fitting from the distribution of amplitudes in the 2412 wavenumber domain. 2413 both float and string values accepted. 2414 given a float value, the unit is set to sigma. 2415 for string values, allowed formats include: 2416 'xsigma' or 'x' (= x-sigma level. e.g., '3sigma'), or 2417 'topx' (= the x strongest ones, e.g. 'top5'). 2418 default is 3.0 (unit: sigma). 2419 addwn: the additional wave numbers to be used for fitting. 2420 list or integer value is accepted to specify every 2421 wave numbers. also string value can be used in case 2422 you need to specify wave numbers in a certain range, 2423 e.g., 'a-b' (= a, a+1, a+2, ..., b-1, b), 2424 '<a' (= 0,1,...,a-2,a-1), 2425 '>=a' (= a, a+1, ... up to the maximum wave 2426 number corresponding to the Nyquist 2427 frequency for the case of FFT). 2428 default is []. 2429 rejwn: the wave numbers NOT to be used for fitting. 2430 can be set just as addwn but has higher priority: 2431 wave numbers which are specified both in addwn 2432 and rejwn will NOT be used. default is []. 2250 2433 clipthresh: Clipping threshold. (default is 3.0, unit: sigma) 2251 2434 clipniter: maximum number of iteration of 'clipthresh'-sigma clipping (default is 0) … … 2283 2466 2284 2467 Example: 2285 bscan = scan.auto_sinusoid_baseline( nwave=10, insitu=False)2468 bscan = scan.auto_sinusoid_baseline(addwn='<=10', insitu=False) 2286 2469 2287 2470 Note: … … 2290 2473 """ 2291 2474 2292 varlist = vars() 2293 2294 if insitu is None: insitu = rcParams['insitu'] 2295 if insitu: 2296 workscan = self 2297 else: 2298 workscan = self.copy() 2299 2300 nchan = workscan.nchan() 2301 2302 if mask is None: mask = [True for i in xrange(nchan)] 2303 if nwave is None: nwave = 3 2304 if maxwavelength is None: maxwavelength = 1.0 2305 if clipthresh is None: clipthresh = 3.0 2306 if clipniter is None: clipniter = 0 2307 if edge is None: edge = (0,0) 2308 if threshold is None: threshold = 3 2309 if chan_avg_limit is None: chan_avg_limit = 1 2310 if plot is None: plot = False 2311 if getresidual is None: getresidual = True 2312 if outlog is None: outlog = False 2313 if blfile is None: blfile = "" 2314 2315 outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile))) 2316 2317 from asap.asaplinefind import linefinder 2318 from asap import _is_sequence_or_number as _is_valid 2319 2320 if isinstance(nwave, int): 2321 in_nwave = nwave 2322 nwave = [] 2323 for i in xrange(in_nwave+1): nwave.append(i) 2324 2325 if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ] 2326 individualedge = False; 2327 if len(edge) > 1: individualedge = isinstance(edge[0], list) or isinstance(edge[0], tuple) 2328 2329 if individualedge: 2330 for edgepar in edge: 2331 if not _is_valid(edgepar, int): 2332 raise ValueError, "Each element of the 'edge' tuple has \ 2333 to be a pair of integers or an integer." 2334 else: 2335 if not _is_valid(edge, int): 2336 raise ValueError, "Parameter 'edge' has to be an integer or a \ 2337 pair of integers specified as a tuple. \ 2338 Nested tuples are allowed \ 2339 to make individual selection for different IFs." 2340 2341 if len(edge) > 1: 2342 curedge = edge 2475 try: 2476 varlist = vars() 2477 2478 if insitu is None: insitu = rcParams['insitu'] 2479 if insitu: 2480 workscan = self 2343 2481 else: 2344 curedge = edge + edge 2345 2346 try: 2482 workscan = self.copy() 2483 2484 if mask is None: mask = [True for i in xrange(workscan.nchan())] 2485 if applyfft is None: applyfft = True 2486 if fftmethod is None: fftmethod = 'fft' 2487 if fftthresh is None: fftthresh = 3.0 2488 if addwn is None: addwn = [] 2489 if rejwn is None: rejwn = [] 2490 if clipthresh is None: clipthresh = 3.0 2491 if clipniter is None: clipniter = 0 2492 if edge is None: edge = (0,0) 2493 if threshold is None: threshold = 3 2494 if chan_avg_limit is None: chan_avg_limit = 1 2495 if plot is None: plot = False 2496 if getresidual is None: getresidual = True 2497 if outlog is None: outlog = False 2498 if blfile is None: blfile = '' 2499 2347 2500 #CURRENTLY, PLOT=true is UNAVAILABLE UNTIL sinusoidal fitting is implemented as a fitter method. 2348 if individualedge: 2349 curedge = [] 2350 for i in xrange(len(edge)): 2351 curedge += edge[i] 2352 2353 workscan._auto_sinusoid_baseline(mask, nwave, maxwavelength, clipthresh, clipniter, curedge, threshold, chan_avg_limit, getresidual, outlog, blfile) 2354 2501 workscan._auto_sinusoid_baseline(mask, applyfft, fftmethod.lower(), str(fftthresh).lower(), workscan._parse_wn(addwn), workscan._parse_wn(rejwn), clipthresh, clipniter, normalise_edge_param(edge), threshold, chan_avg_limit, getresidual, outlog, blfile) 2355 2502 workscan._add_history("auto_sinusoid_baseline", varlist) 2356 2503 … … 2361 2508 2362 2509 except RuntimeError, e: 2363 msg = "The fit failed, possibly because it didn't converge." 2364 if rcParams["verbose"]: 2365 asaplog.push(str(e)) 2366 asaplog.push(str(msg)) 2367 return 2368 else: 2369 raise RuntimeError(str(e)+'\n'+msg) 2370 2510 raise_fitting_failure_exception(e) 2371 2511 2372 2512 @asaplog_post_dec … … 2405 2545 """ 2406 2546 2407 varlist = vars()2408 2409 if insitu is None: insitu = rcParams["insitu"]2410 if insitu:2411 workscan = self2412 else:2413 workscan = self.copy()2414 2415 nchan = workscan.nchan()2416 2417 if mask is None: mask = [True for i in xrange(nchan)]2418 if npiece is None: npiece = 22419 if clipthresh is None: clipthresh = 3.02420 if clipniter is None: clipniter = 02421 if plot is None: plot = False2422 if getresidual is None: getresidual = True2423 if outlog is None: outlog = False2424 if blfile is None: blfile = ""2425 2426 outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))2427 2428 2547 try: 2548 varlist = vars() 2549 2550 if insitu is None: insitu = rcParams['insitu'] 2551 if insitu: 2552 workscan = self 2553 else: 2554 workscan = self.copy() 2555 2556 if mask is None: mask = [True for i in xrange(workscan.nchan())] 2557 if npiece is None: npiece = 2 2558 if clipthresh is None: clipthresh = 3.0 2559 if clipniter is None: clipniter = 0 2560 if plot is None: plot = False 2561 if getresidual is None: getresidual = True 2562 if outlog is None: outlog = False 2563 if blfile is None: blfile = '' 2564 2429 2565 #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method. 2430 2566 workscan._cspline_baseline(mask, npiece, clipthresh, clipniter, getresidual, outlog, blfile) 2431 2432 2567 workscan._add_history("cspline_baseline", varlist) 2433 2568 … … 2438 2573 2439 2574 except RuntimeError, e: 2440 msg = "The fit failed, possibly because it didn't converge." 2441 if rcParams["verbose"]: 2442 asaplog.push(str(e)) 2443 asaplog.push(str(msg)) 2444 return 2445 else: 2446 raise RuntimeError(str(e)+'\n'+msg) 2447 2448 2575 raise_fitting_failure_exception(e) 2576 2577 @asaplog_post_dec 2449 2578 def auto_cspline_baseline(self, insitu=None, mask=None, npiece=None, clipthresh=None, 2450 2579 clipniter=None, edge=None, threshold=None, … … 2505 2634 """ 2506 2635 2507 varlist = vars() 2508 2509 if insitu is None: insitu = rcParams['insitu'] 2510 if insitu: 2511 workscan = self 2512 else: 2513 workscan = self.copy() 2514 2515 nchan = workscan.nchan() 2516 2517 if mask is None: mask = [True for i in xrange(nchan)] 2518 if npiece is None: npiece = 2 2519 if clipthresh is None: clipthresh = 3.0 2520 if clipniter is None: clipniter = 0 2521 if edge is None: edge = (0, 0) 2522 if threshold is None: threshold = 3 2523 if chan_avg_limit is None: chan_avg_limit = 1 2524 if plot is None: plot = False 2525 if getresidual is None: getresidual = True 2526 if outlog is None: outlog = False 2527 if blfile is None: blfile = "" 2528 2529 outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile))) 2530 2531 from asap.asaplinefind import linefinder 2532 from asap import _is_sequence_or_number as _is_valid 2533 2534 if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ] 2535 individualedge = False; 2536 if len(edge) > 1: individualedge = isinstance(edge[0], list) or isinstance(edge[0], tuple) 2537 2538 if individualedge: 2539 for edgepar in edge: 2540 if not _is_valid(edgepar, int): 2541 raise ValueError, "Each element of the 'edge' tuple has \ 2542 to be a pair of integers or an integer." 2543 else: 2544 if not _is_valid(edge, int): 2545 raise ValueError, "Parameter 'edge' has to be an integer or a \ 2546 pair of integers specified as a tuple. \ 2547 Nested tuples are allowed \ 2548 to make individual selection for different IFs." 2549 2550 if len(edge) > 1: 2551 curedge = edge 2636 try: 2637 varlist = vars() 2638 2639 if insitu is None: insitu = rcParams['insitu'] 2640 if insitu: 2641 workscan = self 2552 2642 else: 2553 curedge = edge + edge 2554 2555 try: 2556 #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method. 2557 if individualedge: 2558 curedge = [] 2559 for i in xrange(len(edge)): 2560 curedge += edge[i] 2561 2562 workscan._auto_cspline_baseline(mask, npiece, clipthresh, clipniter, curedge, threshold, chan_avg_limit, getresidual, outlog, blfile) 2563 2643 workscan = self.copy() 2644 2645 if mask is None: mask = [True for i in xrange(workscan.nchan())] 2646 if npiece is None: npiece = 2 2647 if clipthresh is None: clipthresh = 3.0 2648 if clipniter is None: clipniter = 0 2649 if edge is None: edge = (0, 0) 2650 if threshold is None: threshold = 3 2651 if chan_avg_limit is None: chan_avg_limit = 1 2652 if plot is None: plot = False 2653 if getresidual is None: getresidual = True 2654 if outlog is None: outlog = False 2655 if blfile is None: blfile = '' 2656 2657 #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method. 2658 workscan._auto_cspline_baseline(mask, npiece, clipthresh, clipniter, normalise_edge_param(edge), threshold, chan_avg_limit, getresidual, outlog, blfile) 2564 2659 workscan._add_history("auto_cspline_baseline", varlist) 2565 2660 … … 2570 2665 2571 2666 except RuntimeError, e: 2572 msg = "The fit failed, possibly because it didn't converge." 2573 if rcParams["verbose"]: 2574 asaplog.push(str(e)) 2575 asaplog.push(str(msg)) 2576 return 2577 else: 2578 raise RuntimeError(str(e)+'\n'+msg) 2579 2667 raise_fitting_failure_exception(e) 2580 2668 2581 2669 @asaplog_post_dec … … 2606 2694 """ 2607 2695 2608 varlist = vars() 2696 try: 2697 varlist = vars() 2609 2698 2610 if insitu is None: insitu = rcParams["insitu"] 2611 if insitu: 2612 workscan = self 2613 else: 2614 workscan = self.copy() 2615 2616 nchan = workscan.nchan() 2617 2618 if mask is None: mask = [True for i in xrange(nchan)] 2619 if order is None: order = 0 2620 if plot is None: plot = False 2621 if getresidual is None: getresidual = True 2622 if outlog is None: outlog = False 2623 if blfile is None: blfile = "" 2624 2625 outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile))) 2626 2627 try: 2628 rows = xrange(workscan.nrow()) 2629 2630 #if len(rows) > 0: workscan._init_blinfo() 2699 if insitu is None: insitu = rcParams["insitu"] 2700 if insitu: 2701 workscan = self 2702 else: 2703 workscan = self.copy() 2704 2705 if mask is None: mask = [True for i in xrange(workscan.nchan())] 2706 if order is None: order = 0 2707 if plot is None: plot = False 2708 if getresidual is None: getresidual = True 2709 if outlog is None: outlog = False 2710 if blfile is None: blfile = "" 2631 2711 2632 2712 if plot: 2713 outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile))) 2633 2714 if outblfile: blf = open(blfile, "a") 2634 2715 2635 2716 f = fitter() 2636 2717 f.set_function(lpoly=order) 2718 2719 rows = xrange(workscan.nrow()) 2720 #if len(rows) > 0: workscan._init_blinfo() 2721 2637 2722 for r in rows: 2638 2723 f.x = workscan._getabcissa(r) … … 2673 2758 2674 2759 except RuntimeError, e: 2675 msg = "The fit failed, possibly because it didn't converge." 2676 if rcParams["verbose"]: 2677 asaplog.push(str(e)) 2678 asaplog.push(str(msg)) 2679 return 2680 else: 2681 raise RuntimeError(str(e)+'\n'+msg) 2682 2683 2760 raise_fitting_failure_exception(e) 2761 2762 @asaplog_post_dec 2684 2763 def auto_poly_baseline(self, insitu=None, mask=None, order=None, edge=None, threshold=None, 2685 2764 chan_avg_limit=None, plot=None, getresidual=None, outlog=None, blfile=None): … … 2731 2810 """ 2732 2811 2733 varlist = vars() 2734 2735 if insitu is None: insitu = rcParams['insitu'] 2736 if insitu: 2737 workscan = self 2738 else: 2739 workscan = self.copy() 2740 2741 nchan = workscan.nchan() 2742 2743 if mask is None: mask = [True for i in xrange(nchan)] 2744 if order is None: order = 0 2745 if edge is None: edge = (0, 0) 2746 if threshold is None: threshold = 3 2747 if chan_avg_limit is None: chan_avg_limit = 1 2748 if plot is None: plot = False 2749 if getresidual is None: getresidual = True 2750 if outlog is None: outlog = False 2751 if blfile is None: blfile = "" 2752 2753 outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile))) 2754 2755 from asap.asaplinefind import linefinder 2756 from asap import _is_sequence_or_number as _is_valid 2757 2758 if not (isinstance(edge, list) or isinstance(edge, tuple)): edge = [ edge ] 2759 individualedge = False; 2760 if len(edge) > 1: individualedge = isinstance(edge[0], list) or isinstance(edge[0], tuple) 2761 2762 if individualedge: 2763 for edgepar in edge: 2764 if not _is_valid(edgepar, int): 2765 raise ValueError, "Each element of the 'edge' tuple has \ 2766 to be a pair of integers or an integer." 2767 else: 2768 if not _is_valid(edge, int): 2769 raise ValueError, "Parameter 'edge' has to be an integer or a \ 2770 pair of integers specified as a tuple. \ 2771 Nested tuples are allowed \ 2772 to make individual selection for different IFs." 2773 2774 if len(edge) > 1: 2775 curedge = edge 2812 try: 2813 varlist = vars() 2814 2815 if insitu is None: insitu = rcParams['insitu'] 2816 if insitu: 2817 workscan = self 2776 2818 else: 2777 curedge = edge + edge 2778 2779 try: 2780 rows = xrange(workscan.nrow()) 2781 2782 #if len(rows) > 0: workscan._init_blinfo() 2819 workscan = self.copy() 2820 2821 if mask is None: mask = [True for i in xrange(workscan.nchan())] 2822 if order is None: order = 0 2823 if edge is None: edge = (0, 0) 2824 if threshold is None: threshold = 3 2825 if chan_avg_limit is None: chan_avg_limit = 1 2826 if plot is None: plot = False 2827 if getresidual is None: getresidual = True 2828 if outlog is None: outlog = False 2829 if blfile is None: blfile = '' 2830 2831 edge = normalise_edge_param(edge) 2783 2832 2784 2833 if plot: 2834 outblfile = (blfile != "") and os.path.exists(os.path.expanduser(os.path.expandvars(blfile))) 2785 2835 if outblfile: blf = open(blfile, "a") 2786 2836 2837 from asap.asaplinefind import linefinder 2787 2838 fl = linefinder() 2788 2839 fl.set_options(threshold=threshold,avg_limit=chan_avg_limit) 2789 2840 fl.set_scan(workscan) 2841 2790 2842 f = fitter() 2791 2843 f.set_function(lpoly=order) 2792 2844 2845 rows = xrange(workscan.nrow()) 2846 #if len(rows) > 0: workscan._init_blinfo() 2847 2793 2848 for r in rows: 2794 if individualedge: 2795 if len(edge) <= workscan.getif(r): 2796 raise RuntimeError, "Number of edge elements appear to " \ 2797 "be less than the number of IFs" 2798 else: 2799 curedge = edge[workscan.getif(r)] 2800 2801 fl.find_lines(r, mask_and(mask, workscan._getmask(r)), curedge) # (CAS-1434) 2849 idx = 2*workscan.getif(r) 2850 fl.find_lines(r, mask_and(mask, workscan._getmask(r)), edge[idx:idx+2]) # (CAS-1434) 2802 2851 2803 2852 f.x = workscan._getabcissa(r) … … 2827 2876 2828 2877 if outblfile: blf.close() 2829 2830 2878 else: 2831 if individualedge: 2832 curedge = [] 2833 for i in xrange(len(edge)): 2834 curedge += edge[i] 2835 2836 workscan._auto_poly_baseline(mask, order, curedge, threshold, chan_avg_limit, getresidual, outlog, blfile) 2879 workscan._auto_poly_baseline(mask, order, edge, threshold, chan_avg_limit, getresidual, outlog, blfile) 2837 2880 2838 2881 workscan._add_history("auto_poly_baseline", varlist) … … 2844 2887 2845 2888 except RuntimeError, e: 2846 msg = "The fit failed, possibly because it didn't converge." 2847 if rcParams["verbose"]: 2848 asaplog.push(str(e)) 2849 asaplog.push(str(msg)) 2850 return 2851 else: 2852 raise RuntimeError(str(e)+'\n'+msg) 2853 2854 2855 ### OBSOLETE ################################################################## 2856 @asaplog_post_dec 2857 def old_poly_baseline(self, mask=None, order=0, plot=False, uselin=False, insitu=None, rows=None): 2858 """ 2859 Return a scan which has been baselined (all rows) by a polynomial. 2860 2861 Parameters: 2862 2863 mask: an optional mask 2864 2865 order: the order of the polynomial (default is 0) 2866 2867 plot: plot the fit and the residual. In this each 2868 indivual fit has to be approved, by typing 'y' 2869 or 'n' 2870 2871 uselin: use linear polynomial fit 2872 2873 insitu: if False a new scantable is returned. 2874 Otherwise, the scaling is done in-situ 2875 The default is taken from .asaprc (False) 2876 2877 rows: row numbers of spectra to be processed. 2878 (default is None: for all rows) 2879 2880 Example: 2881 # return a scan baselined by a third order polynomial, 2882 # not using a mask 2883 bscan = scan.poly_baseline(order=3) 2884 2885 """ 2886 if insitu is None: insitu = rcParams['insitu'] 2887 if not insitu: 2888 workscan = self.copy() 2889 else: 2890 workscan = self 2891 varlist = vars() 2892 if mask is None: 2893 mask = [True for i in xrange(self.nchan())] 2894 2895 try: 2896 f = fitter() 2897 if uselin: 2898 f.set_function(lpoly=order) 2899 else: 2900 f.set_function(poly=order) 2901 2902 if rows == None: 2903 rows = xrange(workscan.nrow()) 2904 elif isinstance(rows, int): 2905 rows = [ rows ] 2906 2907 if len(rows) > 0: 2908 self.blpars = [] 2909 self.masklists = [] 2910 self.actualmask = [] 2911 2912 for r in rows: 2913 f.x = workscan._getabcissa(r) 2914 f.y = workscan._getspectrum(r) 2915 f.mask = mask_and(mask, workscan._getmask(r)) # (CAS-1434) 2916 f.data = None 2917 f.fit() 2918 if plot: 2919 f.plot(residual=True) 2920 x = raw_input("Accept fit ( [y]/n ): ") 2921 if x.upper() == 'N': 2922 self.blpars.append(None) 2923 self.masklists.append(None) 2924 self.actualmask.append(None) 2925 continue 2926 workscan._setspectrum(f.fitter.getresidual(), r) 2927 self.blpars.append(f.get_parameters()) 2928 self.masklists.append(workscan.get_masklist(f.mask, row=r, silent=True)) 2929 self.actualmask.append(f.mask) 2930 2931 if plot: 2932 f._p.unmap() 2933 f._p = None 2934 workscan._add_history("poly_baseline", varlist) 2935 if insitu: 2936 self._assign(workscan) 2937 else: 2938 return workscan 2939 except RuntimeError: 2940 msg = "The fit failed, possibly because it didn't converge." 2941 raise RuntimeError(msg) 2889 raise_fitting_failure_exception(e) 2942 2890 2943 2891 def _init_blinfo(self):
Note:
See TracChangeset
for help on using the changeset viewer.