|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 import os, shutil 00002 import numpy 00003 import numpy.fft as FFT 00004 import math 00005 00006 from asap.scantable import scantable 00007 from asap.parameters import rcParams 00008 from asap.logging import asaplog, asaplog_post_dec 00009 from asap.selector import selector 00010 from asap.asapgrid import asapgrid2 00011 #from asap._asap import sidebandsep 00012 00013 class sbseparator: 00014 """ 00015 The sbseparator class is defined to separate SIGNAL and IMAGE 00016 sideband spectra observed by frequency-switching technique. 00017 It also helps supressing emmission of IMAGE sideband. 00018 *** WARNING *** THIS MODULE IS EXPERIMENTAL 00019 Known issues: 00020 - Frequencies of IMAGE sideband cannot be reconstructed from 00021 information in scantable in sideband sparation. Frequency 00022 setting of SIGNAL sideband is stored in output table for now. 00023 - Flag information (stored in FLAGTRA) is ignored. 00024 00025 Example: 00026 # Create sideband separator instance whith 3 input data 00027 sbsep = sbseparator(['test1.asap', 'test2.asap', 'test3.asap']) 00028 # Set reference IFNO and tolerance to select data 00029 sbsep.set_frequency(5, 30, frame='TOPO') 00030 # Set direction tolerance to select data in unit of radian 00031 sbsep.set_dirtol(1.e-5) 00032 # Set rejection limit of solution 00033 sbsep.set_limit(0.2) 00034 # Solve image sideband as well 00035 sbsep.set_both(True) 00036 # Invoke sideband separation 00037 sbsep.separate('testout.asap', overwrite = True) 00038 """ 00039 def __init__(self, infiles): 00040 self.intables = None 00041 self.signalShift = [] 00042 self.imageShift = [] 00043 self.dsbmode = True 00044 self.getboth = False 00045 self.rejlimit = 0.2 00046 self.baseif = -1 00047 self.freqtol = 10. 00048 self.freqframe = "" 00049 self.solveother = False 00050 self.dirtol = [1.e-5, 1.e-5] # direction tolerance in rad (2 arcsec) 00051 00052 self.tables = [] 00053 self.nshift = -1 00054 self.nchan = -1 00055 00056 self.set_data(infiles) 00057 00058 #self.separator = sidebandsep() 00059 00060 @asaplog_post_dec 00061 def set_data(self, infiles): 00062 """ 00063 Set data to be processed. 00064 00065 infiles : a list of filenames or scantables 00066 """ 00067 if not (type(infiles) in (list, tuple, numpy.ndarray)): 00068 infiles = [infiles] 00069 if isinstance(infiles[0], scantable): 00070 # a list of scantable 00071 for stab in infiles: 00072 if not isinstance(stab, scantable): 00073 asaplog.post() 00074 raise TypeError, "Input data is not a list of scantables." 00075 00076 #self.separator._setdata(infiles) 00077 self._reset_data() 00078 self.intables = infiles 00079 else: 00080 # a list of filenames 00081 for name in infiles: 00082 if not os.path.exists(name): 00083 asaplog.post() 00084 raise ValueError, "Could not find input file '%s'" % name 00085 00086 #self.separator._setdataname(infiles) 00087 self._reset_data() 00088 self.intables = infiles 00089 00090 asaplog.push("%d files are set to process" % len(self.intables)) 00091 00092 00093 def _reset_data(self): 00094 del self.intables 00095 self.intables = None 00096 self.signalShift = [] 00097 #self.imageShift = [] 00098 self.tables = [] 00099 self.nshift = -1 00100 self.nchan = -1 00101 00102 @asaplog_post_dec 00103 def set_frequency(self, baseif, freqtol, frame=""): 00104 """ 00105 Set IFNO and frequency tolerance to select data to process. 00106 00107 Parameters: 00108 - reference IFNO to process in the first table in the list 00109 - frequency tolerance from reference IF to select data 00110 frame : frequency frame to select IF 00111 """ 00112 self._reset_if() 00113 self.baseif = baseif 00114 if isinstance(freqtol,dict) and freqtol["unit"] == "Hz": 00115 if freqtol['value'] > 0.: 00116 self.freqtol = freqtol 00117 else: 00118 asaplog.post() 00119 asaplog.push("Frequency tolerance should be positive value.") 00120 asaplog.post("ERROR") 00121 return 00122 else: 00123 # torelance in channel unit 00124 if freqtol > 0: 00125 self.freqtol = float(freqtol) 00126 else: 00127 asaplog.post() 00128 asaplog.push("Frequency tolerance should be positive value.") 00129 asaplog.post("ERROR") 00130 return 00131 self.freqframe = frame 00132 00133 def _reset_if(self): 00134 self.baseif = -1 00135 self.freqtol = 0 00136 self.freqframe = "" 00137 self.signalShift = [] 00138 #self.imageShift = [] 00139 self.tables = [] 00140 self.nshift = 0 00141 self.nchan = -1 00142 00143 @asaplog_post_dec 00144 def set_dirtol(self, dirtol=[1.e-5,1.e-5]): 00145 """ 00146 Set tolerance of direction to select data 00147 """ 00148 # direction tolerance in rad 00149 if not (type(dirtol) in [list, tuple, numpy.ndarray]): 00150 dirtol = [dirtol, dirtol] 00151 if len(dirtol) == 1: 00152 dirtol = [dirtol[0], dirtol[0]] 00153 if len(dirtol) > 1: 00154 self.dirtol = dirtol[0:2] 00155 else: 00156 asaplog.post() 00157 asaplog.push("Invalid direction tolerance. Should be a list of float in unit radian") 00158 asaplog.post("ERROR") 00159 return 00160 asaplog.post("Set direction tolerance [%f, %f] (rad)" % \ 00161 (self.dirtol[0], self.dirtol[1])) 00162 00163 @asaplog_post_dec 00164 def set_shift(self, mode="DSB", imageshift=None): 00165 """ 00166 Set shift mode and channel shift of image band. 00167 00168 mode : shift mode ['DSB'|'SSB'] 00169 When mode='DSB', imageshift is assumed to be equal 00170 to the shift of signal sideband but in opposite direction. 00171 imageshift : a list of number of channel shift in image band of 00172 each scantable. valid only mode='SSB' 00173 """ 00174 if mode.upper().startswith("S"): 00175 if not imageshift: 00176 raise ValueError, "Need to set shift value of image sideband" 00177 self.dsbmode = False 00178 self.imageShift = imageshift 00179 asaplog.push("Image sideband shift is set manually: %s" % str(self.imageShift)) 00180 else: 00181 # DSB mode 00182 self.dsbmode = True 00183 self.imageShift = [] 00184 00185 @asaplog_post_dec 00186 def set_both(self, flag=False): 00187 """ 00188 Resolve both image and signal sideband when True. 00189 """ 00190 self.getboth = flag 00191 if self.getboth: 00192 asaplog.push("Both signal and image sidebands are solved and output as separate tables.") 00193 else: 00194 asaplog.push("Only signal sideband is solved and output as an table.") 00195 00196 @asaplog_post_dec 00197 def set_limit(self, threshold=0.2): 00198 """ 00199 Set rejection limit of solution. 00200 """ 00201 #self.separator._setlimit(abs(threshold)) 00202 self.rejlimit = threshold 00203 asaplog.push("The threshold of rejection is set to %f" % self.rejlimit) 00204 00205 00206 @asaplog_post_dec 00207 def set_solve_other(self, flag=False): 00208 """ 00209 Calculate spectra by subtracting the solution of the other sideband 00210 when True. 00211 """ 00212 self.solveother = flag 00213 if flag: 00214 asaplog.push("Expert mode: solution are obtained by subtraction of the other sideband.") 00215 00216 00217 @asaplog_post_dec 00218 def separate(self, outname="", overwrite=False): 00219 """ 00220 Invoke sideband separation. 00221 00222 outname : a name of output scantable 00223 overwrite : overwrite existing table 00224 """ 00225 # List up valid scantables and IFNOs to convolve. 00226 #self.separator._separate() 00227 self._setup_shift() 00228 #self._preprocess_tables() 00229 00230 nshift = len(self.tables) 00231 signaltab = self._grid_outtable(self.tables[0].copy()) 00232 if self.getboth: 00233 imagetab = signaltab.copy() 00234 00235 rejrow = [] 00236 for irow in xrange(signaltab.nrow()): 00237 currpol = signaltab.getpol(irow) 00238 currbeam = signaltab.getbeam(irow) 00239 currdir = signaltab.get_directionval(irow) 00240 spec_array, tabidx = self._get_specarray(polid=currpol,\ 00241 beamid=currbeam,\ 00242 dir=currdir) 00243 #if not spec_array: 00244 if len(tabidx) == 0: 00245 asaplog.post() 00246 asaplog.push("skipping row %d" % irow) 00247 rejrow.append(irow) 00248 continue 00249 signal = self._solve_signal(spec_array, tabidx) 00250 signaltab.set_spectrum(signal, irow) 00251 if self.getboth: 00252 image = self._solve_image(spec_array, tabidx) 00253 imagetab.set_spectrum(image, irow) 00254 00255 # TODO: Need to remove rejrow form scantables here 00256 signaltab.flag_row(rejrow) 00257 if self.getboth: 00258 imagetab.flag_row(rejrow) 00259 00260 if outname == "": 00261 outname = "sbsepareted.asap" 00262 signalname = outname + ".signalband" 00263 if os.path.exists(signalname): 00264 if not overwrite: 00265 raise Exception, "Output file '%s' exists." % signalname 00266 else: 00267 shutil.rmtree(signalname) 00268 signaltab.save(signalname) 00269 if self.getboth: 00270 # Warnings 00271 asaplog.post() 00272 asaplog.push("Saving IMAGE sideband.") 00273 asaplog.push("Note, frequency information of IMAGE sideband cannot be properly filled so far. (future development)") 00274 asaplog.push("Storing frequency setting of SIGNAL sideband in FREQUENCIES table for now.") 00275 asaplog.post("WARN") 00276 00277 imagename = outname + ".imageband" 00278 if os.path.exists(imagename): 00279 if not overwrite: 00280 raise Exception, "Output file '%s' exists." % imagename 00281 else: 00282 shutil.rmtree(imagename) 00283 imagetab.save(imagename) 00284 00285 00286 def _solve_signal(self, data, tabidx=None): 00287 if not tabidx: 00288 tabidx = range(len(data)) 00289 00290 tempshift = [] 00291 dshift = [] 00292 if self.solveother: 00293 for idx in tabidx: 00294 tempshift.append(-self.imageShift[idx]) 00295 dshift.append(self.signalShift[idx] - self.imageShift[idx]) 00296 else: 00297 for idx in tabidx: 00298 tempshift.append(-self.signalShift[idx]) 00299 dshift.append(self.imageShift[idx] - self.signalShift[idx]) 00300 00301 shiftdata = numpy.zeros(data.shape, numpy.float) 00302 for i in range(len(data)): 00303 shiftdata[i] = self._shiftSpectrum(data[i], tempshift[i]) 00304 ifftdata = self._Deconvolution(shiftdata, dshift, self.rejlimit) 00305 result_image = self._combineResult(ifftdata) 00306 if not self.solveother: 00307 return result_image 00308 result_signal = self._subtractOtherSide(shiftdata, dshift, result_image) 00309 return result_signal 00310 00311 00312 def _solve_image(self, data, tabidx=None): 00313 if not tabidx: 00314 tabidx = range(len(data)) 00315 00316 tempshift = [] 00317 dshift = [] 00318 if self.solveother: 00319 for idx in tabidx: 00320 tempshift.append(-self.signalShift[idx]) 00321 dshift.append(self.imageShift[idx] - self.signalShift[idx]) 00322 else: 00323 for idx in tabidx: 00324 tempshift.append(-self.imageShift[idx]) 00325 dshift.append(self.signalShift[idx] - self.imageShift[idx]) 00326 00327 shiftdata = numpy.zeros(data.shape, numpy.float) 00328 for i in range(len(data)): 00329 shiftdata[i] = self._shiftSpectrum(data[i], tempshift[i]) 00330 ifftdata = self._Deconvolution(shiftdata, dshift, self.rejlimit) 00331 result_image = self._combineResult(ifftdata) 00332 if not self.solveother: 00333 return result_image 00334 result_signal = self._subtractOtherSide(shiftdata, dshift, result_image) 00335 return result_signal 00336 00337 @asaplog_post_dec 00338 def _grid_outtable(self, table): 00339 # Generate gridded table for output (Just to get rows) 00340 gridder = asapgrid2(table) 00341 gridder.setIF(self.baseif) 00342 00343 cellx = str(self.dirtol[0])+"rad" 00344 celly = str(self.dirtol[1])+"rad" 00345 dirarr = numpy.array(table.get_directionval()).transpose() 00346 mapx = dirarr[0].max() - dirarr[0].min() 00347 mapy = dirarr[1].max() - dirarr[1].min() 00348 nx = max(1, numpy.ceil(mapx/self.dirtol[0])) 00349 ny = max(1, numpy.ceil(mapy/self.dirtol[0])) 00350 00351 asaplog.push("Regrid output scantable with cell = [%s, %s]" % \ 00352 (cellx, celly)) 00353 gridder.defineImage(nx=nx, ny=ny, cellx=cellx, celly=celly) 00354 gridder.setFunc(func='box', width=1) 00355 gridder.setWeight(weightType='uniform') 00356 gridder.grid() 00357 return gridder.getResult() 00358 00359 00360 @asaplog_post_dec 00361 def _get_specarray(self, polid=None, beamid=None, dir=None): 00362 ntable = len(self.tables) 00363 spec_array = numpy.zeros((ntable, self.nchan), numpy.float) 00364 nspec = 0 00365 asaplog.push("Start data selection by POL=%d, BEAM=%d, direction=[%f, %f]" % (polid, beamid, dir[0], dir[1])) 00366 tabidx = [] 00367 for itab in range(ntable): 00368 tab = self.tables[itab] 00369 # Select rows by POLNO and BEAMNO 00370 try: 00371 tab.set_selection(pols=[polid], beams=[beamid]) 00372 if tab.nrow() > 0: tabidx.append(itab) 00373 except: # no selection 00374 asaplog.post() 00375 asaplog.push("table %d - No spectrum ....skipping the table" % (itab)) 00376 asaplog.post("WARN") 00377 continue 00378 00379 # Select rows by direction 00380 spec = numpy.zeros(self.nchan, numpy.float) 00381 selrow = [] 00382 for irow in range(tab.nrow()): 00383 currdir = tab.get_directionval(irow) 00384 if (abs(currdir[0] - dir[0]) > self.dirtol[0]) or \ 00385 (abs(currdir[1] - dir[1]) > self.dirtol[1]): 00386 continue 00387 selrow.append(irow) 00388 if len(selrow) == 0: 00389 asaplog.post() 00390 asaplog.push("table %d - No spectrum ....skipping the table" % (itab)) 00391 asaplog.post("WARN") 00392 continue 00393 else: 00394 seltab = tab.copy() 00395 seltab.set_selection(selector(rows=selrow)) 00396 00397 if tab.nrow() > 1: 00398 asaplog.push("table %d - More than a spectrum selected. averaging rows..." % (itab)) 00399 tab = seltab.average_time(scanav=False, weight="tintsys") 00400 else: 00401 tab = seltab 00402 00403 spec_array[nspec] = tab._getspectrum() 00404 nspec += 1 00405 00406 if nspec != ntable: 00407 asaplog.post() 00408 #asaplog.push("Some tables has no spectrum with POL=%d BEAM=%d. averaging rows..." % (polid, beamid)) 00409 asaplog.push("Could not find corresponding rows in some tables.") 00410 asaplog.push("Number of spectra selected = %d (table: %d)" % (nspec, ntable)) 00411 if nspec < 2: 00412 asaplog.push("At least 2 spectra are necessary for convolution") 00413 asaplog.post("ERROR") 00414 return False, tabidx 00415 00416 return spec_array[0:nspec], tabidx 00417 00418 00419 @asaplog_post_dec 00420 def _setup_shift(self): 00421 ### define self.tables, self.signalShift, and self.imageShift 00422 if not self.intables: 00423 asaplog.post() 00424 raise RuntimeError, "Input data is not defined." 00425 #if self.baseif < 0: 00426 # asaplog.post() 00427 # raise RuntimeError, "Reference IFNO is not defined." 00428 00429 byname = False 00430 #if not self.intables: 00431 if isinstance(self.intables[0], str): 00432 # A list of file name is given 00433 if not os.path.exists(self.intables[0]): 00434 asaplog.post() 00435 raise RuntimeError, "Could not find '%s'" % self.intables[0] 00436 00437 stab = scantable(self.intables[0],average=False) 00438 ntab = len(self.intables) 00439 byname = True 00440 else: 00441 stab = self.intables[0] 00442 ntab = len(self.intables) 00443 00444 if len(stab.getbeamnos()) > 1: 00445 asaplog.post() 00446 asaplog.push("Mult-beam data is not supported by this module.") 00447 asaplog.post("ERROR") 00448 return 00449 00450 valid_ifs = stab.getifnos() 00451 if self.baseif < 0: 00452 self.baseif = valid_ifs[0] 00453 asaplog.post() 00454 asaplog.push("IFNO is not selected. Using the first IF in the first scantable. Reference IFNO = %d" % (self.baseif)) 00455 00456 if not (self.baseif in valid_ifs): 00457 asaplog.post() 00458 errmsg = "IF%d does not exist in the first scantable" % \ 00459 self.baseif 00460 raise RuntimeError, errmsg 00461 00462 asaplog.push("Start selecting tables and IFNOs to solve.") 00463 asaplog.push("Cheching frequency of the reference IF") 00464 unit_org = stab.get_unit() 00465 coord = stab._getcoordinfo() 00466 frame_org = coord[1] 00467 stab.set_unit("Hz") 00468 if len(self.freqframe) > 0: 00469 stab.set_freqframe(self.freqframe) 00470 stab.set_selection(ifs=[self.baseif]) 00471 spx = stab._getabcissa() 00472 stab.set_selection() 00473 basech0 = spx[0] 00474 baseinc = spx[1]-spx[0] 00475 self.nchan = len(spx) 00476 if isinstance(self.freqtol, float): 00477 vftol = abs(baseinc * self.freqtol) 00478 self.freqtol = dict(value=vftol, unit="Hz") 00479 else: 00480 vftol = abs(self.freqtol['value']) 00481 inctol = abs(baseinc/float(self.nchan)) 00482 asaplog.push("Reference frequency setup (Table = 0, IFNO = %d): nchan = %d, chan0 = %f Hz, incr = %f Hz" % (self.baseif, self.nchan, basech0, baseinc)) 00483 asaplog.push("Allowed frequency tolerance = %f Hz ( %f channels)" % (vftol, vftol/baseinc)) 00484 poltype0 = stab.poltype() 00485 00486 self.tables = [] 00487 self.signalShift = [] 00488 if self.dsbmode: 00489 self.imageShift = [] 00490 00491 for itab in range(ntab): 00492 asaplog.push("Table %d:" % itab) 00493 tab_selected = False 00494 if itab > 0: 00495 if byname: 00496 stab = scantable(self.intables[itab],average=False) 00497 self.intables.append(stab) 00498 else: 00499 stab = self.intables[itab] 00500 unit_org = stab.get_unit() 00501 coord = stab._getcoordinfo() 00502 frame_org = coord[1] 00503 stab.set_unit("Hz") 00504 if len(self.freqframe) > 0: 00505 stab.set_freqframe(self.freqframe) 00506 00507 # Check POLTYPE should be equal to the first table. 00508 if stab.poltype() != poltype0: 00509 asaplog.post() 00510 raise Exception, "POLTYPE should be equal to the first table." 00511 # Multiple beam data may not handled properly 00512 if len(stab.getbeamnos()) > 1: 00513 asaplog.post() 00514 asaplog.push("table contains multiple beams. It may not be handled properly.") 00515 asaplog.push("WARN") 00516 00517 for ifno in stab.getifnos(): 00518 stab.set_selection(ifs=[ifno]) 00519 spx = stab._getabcissa() 00520 if (len(spx) != self.nchan) or \ 00521 (abs(spx[0]-basech0) > vftol) or \ 00522 (abs(spx[1]-spx[0]-baseinc) > inctol): 00523 continue 00524 tab_selected = True 00525 seltab = stab.copy() 00526 seltab.set_unit(unit_org) 00527 seltab.set_freqframe(frame_org) 00528 self.tables.append(seltab) 00529 self.signalShift.append((spx[0]-basech0)/baseinc) 00530 if self.dsbmode: 00531 self.imageShift.append(-self.signalShift[-1]) 00532 asaplog.push("- IF%d selected: sideband shift = %16.12e channels" % (ifno, self.signalShift[-1])) 00533 stab.set_selection() 00534 stab.set_unit(unit_org) 00535 stab.set_freqframe(frame_org) 00536 if not tab_selected: 00537 asaplog.post() 00538 asaplog.push("No data selected in Table %d" % itab) 00539 asaplog.post("WARN") 00540 00541 asaplog.push("Total number of IFs selected = %d" % len(self.tables)) 00542 if len(self.tables) < 2: 00543 asaplog.post() 00544 raise RuntimeError, "At least 2 IFs are necessary for convolution!" 00545 00546 if not self.dsbmode and len(self.imageShift) != len(self.signalShift): 00547 asaplog.post() 00548 errmsg = "User defined channel shift of image sideband has %d elements, while selected IFNOs are %d" % (len(self.imageShift), len(self.signalShift)) 00549 raise RuntimeError, errmsg 00550 00551 self.signalShift = numpy.array(self.signalShift) 00552 self.imageShift = numpy.array(self.imageShift) 00553 self.nshift = len(self.tables) 00554 00555 @asaplog_post_dec 00556 def _preprocess_tables(self): 00557 ### temporary method to preprocess data 00558 ### Do time averaging for now. 00559 for itab in range(len(self.tables)): 00560 self.tables[itab] = self.tables[itab].average_time(scanav=False, weight="tintsys") 00561 00562 00563 # def save(self, outfile, outform="ASAP", overwrite=False): 00564 # if not overwrite and os.path.exists(outfile): 00565 # raise RuntimeError, "Output file '%s' already exists" % outfile 00566 # 00567 # #self.separator._save(outfile, outform) 00568 00569 # def done(self): 00570 # self.close() 00571 00572 # def close(self): 00573 # pass 00574 # #del self.separator 00575 00576 00577 00578 ######################################################################## 00579 def _Deconvolution(self, data_array, shift, threshold=0.00000001): 00580 FObs = [] 00581 Reject = 0 00582 nshift, nchan = data_array.shape 00583 nspec = nshift*(nshift-1)/2 00584 ifftObs = numpy.zeros((nspec, nchan), numpy.float) 00585 for i in range(nshift): 00586 F = FFT.fft(data_array[i]) 00587 FObs.append(F) 00588 z = 0 00589 for i in range(nshift): 00590 for j in range(i+1, nshift): 00591 Fobs = (FObs[i]+FObs[j])/2.0 00592 dX = (shift[j]-shift[i])*2.0*math.pi/float(self.nchan) 00593 #print 'dX,i,j=',dX,i,j 00594 for k in range(1,self.nchan): 00595 if math.fabs(math.sin(dX*k)) > threshold: 00596 Fobs[k] += ((FObs[i][k]-FObs[j][k])/2.0/(1.0-math.cos(dX*k))*math.sin(dX*k))*1.0j 00597 else: Reject += 1 00598 ifftObs[z] = FFT.ifft(Fobs) 00599 z += 1 00600 print 'Threshold=%s Reject=%d' % (threshold, Reject) 00601 return ifftObs 00602 00603 def _combineResult(self, ifftObs): 00604 nspec = len(ifftObs) 00605 sum = ifftObs[0] 00606 for i in range(1,nspec): 00607 sum += ifftObs[i] 00608 return(sum/float(nspec)) 00609 00610 def _subtractOtherSide(self, data_array, shift, Data): 00611 sum = numpy.zeros(len(Data), numpy.float) 00612 numSP = len(data_array) 00613 for i in range(numSP): 00614 SPsub = data_array[i] - Data 00615 sum += self._shiftSpectrum(SPsub, -shift[i]) 00616 return(sum/float(numSP)) 00617 00618 def _shiftSpectrum(self, data, Shift): 00619 Out = numpy.zeros(self.nchan, numpy.float) 00620 w2 = Shift % 1 00621 w1 = 1.0 - w2 00622 for i in range(self.nchan): 00623 c1 = int((Shift + i) % self.nchan) 00624 c2 = (c1 + 1) % self.nchan 00625 Out[c1] += data[i] * w1 00626 Out[c2] += data[i] * w2 00627 return Out.copy()
1.8.0
