Changeset 2794
- Timestamp:
- 03/15/13 18:23:04 (12 years ago)
- Location:
- trunk
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/python/sbseparator.py
r2772 r2794 38 38 """ 39 39 def __init__(self, infiles): 40 self.intables = None41 self.signalShift = []42 self.imageShift = []43 self.dsbmode = True44 self.getboth = False45 self.rejlimit = 0.246 self.baseif = -147 self.freqtol = 10.48 self.freqframe = ""49 self.solveother = False50 self.dirtol = [1.e-5, 1.e-5] # direction tolerance in rad (2 arcsec)51 #self.lo1 = 0.52 53 self.tables = []54 self.nshift = -155 self.nchan = -156 57 self.set_data(infiles)58 59 40 self._separator = SBSeparator(infiles) 60 41 61 @asaplog_post_dec 62 def set_data(self, infiles): 63 """ 64 Set data to be processed. 65 66 infiles : a list of filenames or scantables 67 """ 68 if not (type(infiles) in (list, tuple, numpy.ndarray)): 69 infiles = [infiles] 70 if isinstance(infiles[0], scantable): 71 # a list of scantable 72 for stab in infiles: 73 if not isinstance(stab, scantable): 74 asaplog.post() 75 raise TypeError, "Input data is not a list of scantables." 76 77 #self._separator._setdata(infiles) 78 self._reset_data() 79 self.intables = infiles 80 else: 81 # a list of filenames 82 for name in infiles: 83 if not os.path.exists(name): 84 asaplog.post() 85 raise ValueError, "Could not find input file '%s'" % name 86 87 #self._separator._setdataname(infiles) 88 self._reset_data() 89 self.intables = infiles 90 91 asaplog.push("%d files are set to process" % len(self.intables)) 92 93 94 def _reset_data(self): 95 del self.intables 96 self.intables = None 97 self.signalShift = [] 98 #self.imageShift = [] 99 self.tables = [] 100 self.nshift = -1 101 self.nchan = -1 102 103 @asaplog_post_dec 42 ################## temp functions ################# 43 # def rfft(self, invec): 44 # print "Python method cppRfft" 45 # self._separator._cpprfft(invec) 46 # return 47 ################## temp functions ################# 48 104 49 def set_frequency(self, baseif, freqtol, frame=""): 105 50 """ … … 108 53 Parameters: 109 54 - reference IFNO to process in the first table in the list 110 - frequency tolerance from reference IF to select data 55 - frequency tolerance from reference IF to select data (string) 111 56 frame : frequency frame to select IF 112 57 """ 113 self._reset_if() 114 self.baseif = baseif 115 if isinstance(freqtol,dict) and freqtol["unit"] == "Hz": 116 if freqtol['value'] > 0.: 117 self.freqtol = freqtol 58 if type(freqtol) in (float, int): 59 freqtol = str(freqtol) 60 elif isinstance(freqtol, dict): 61 try: 62 freqtol = str(freqtol['value']) + freqtol['unit'] 63 except: 64 raise ValueError, "Invalid frequency tolerance." 65 self._separator.set_freq(baseif, freqtol, frame) 66 67 68 def set_dirtol(self, dirtol=["2arcsec", "2arcsec"]): 69 """ 70 Set tolerance of direction to select data 71 """ 72 if isinstance(dirtol, str): 73 dirtol = [dirtol] 74 75 self._separator.set_dirtol(dirtol) 76 77 78 def set_shift(self, imageshift=[]): 79 """ 80 Set shift mode and channel shift of image band. 81 82 imageshift : a list of number of channels shifted in image 83 side band of each scantable. 84 If the shifts are not set, they are assumed to be 85 equal to those of signal side band, but in opposite 86 direction as usual by LO1 offsetting of DSB receivers. 87 """ 88 if not imageshift: 89 imageshift = [] 90 self._separator.set_shift(imageshift) 91 92 93 @asaplog_post_dec 94 def set_both(self, flag=False): 95 """ 96 Resolve both image and signal sideband when True. 97 """ 98 self._separator.solve_both(flag) 99 if flag: 100 asaplog.push("Both signal and image sidebands will be solved and stored in separate tables.") 101 else: 102 asaplog.push("Only signal sideband will be solved and stored in an table.") 103 104 @asaplog_post_dec 105 def set_limit(self, threshold=0.2): 106 """ 107 Set rejection limit of solution. 108 """ 109 self._separator.set_limit(threshold) 110 111 112 @asaplog_post_dec 113 def set_solve_other(self, flag=False): 114 """ 115 Calculate spectra by subtracting the solution of the other sideband 116 when True. 117 """ 118 self._separator.subtract_other(flag) 119 if flag: 120 asaplog.push("Expert mode: solution are obtained by subtraction of the other sideband.") 121 122 123 def set_lo1(self,lo1, frame="TOPO", reftime=-1, refdir=""): 124 """ 125 Set LO1 frequency to calculate frequency of image sideband. 126 127 lo1 : LO1 frequency 128 frame : the frequency frame of LO1 129 reftime : the reference time used in frequency frame conversion. 130 refdir : the reference direction used in frequency frame conversion. 131 """ 132 self._separator.set_lo1(lo1val, frame, reftime, refdir) 133 134 135 def set_lo1root(self, name): 136 """ 137 Set MS name which stores LO1 frequency of signal side band. 138 It is used to calculate frequency of image sideband. 139 140 name : MS name which contains 'ASDM_SPECTRALWINDOW' and 141 'ASDM_RECEIVER' tables. 142 """ 143 self._separator.set_lo1root(name) 144 145 @asaplog_post_dec 146 def separate(self, outname="", overwrite=False): 147 """ 148 Invoke sideband separation. 149 150 outname : a name of output scantable 151 overwrite : overwrite existing table 152 """ 153 out_default = "sbseparated.asap" 154 if len(outname) == 0: 155 outname = out_default 156 asaplog.post() 157 asaplog.push("The output file name is not specified.") 158 asaplog.push("Using default name '%s'" % outname) 159 asaplog.post("WARN") 160 161 if os.path.exists(outname): 162 if overwrite: 163 asaplog.push("removing the old file '%s'" % outname) 164 shutil.rmtree(outname) 118 165 else: 119 166 asaplog.post() 120 asaplog.push(" Frequency tolerance should be positive value.")167 asaplog.push("Output file '%s' already exists." % outname) 121 168 asaplog.post("ERROR") 122 return 123 else: 124 # torelance in channel unit 125 if freqtol > 0: 126 self.freqtol = float(freqtol) 127 else: 128 asaplog.post() 129 asaplog.push("Frequency tolerance should be positive value.") 130 asaplog.post("ERROR") 131 return 132 self.freqframe = frame 133 134 def _reset_if(self): 135 self.baseif = -1 136 self.freqtol = 0 137 self.freqframe = "" 138 self.signalShift = [] 139 #self.imageShift = [] 140 self.tables = [] 141 self.nshift = 0 142 self.nchan = -1 143 144 @asaplog_post_dec 145 def set_dirtol(self, dirtol=[1.e-5,1.e-5]): 146 """ 147 Set tolerance of direction to select data 148 """ 149 # direction tolerance in rad 150 if not (type(dirtol) in [list, tuple, numpy.ndarray]): 151 dirtol = [dirtol, dirtol] 152 if len(dirtol) == 1: 153 dirtol = [dirtol[0], dirtol[0]] 154 if len(dirtol) > 1: 155 self.dirtol = dirtol[0:2] 156 else: 157 asaplog.post() 158 asaplog.push("Invalid direction tolerance. Should be a list of float in unit radian") 159 asaplog.post("ERROR") 160 return 161 asaplog.post("Set direction tolerance [%f, %f] (rad)" % \ 162 (self.dirtol[0], self.dirtol[1])) 163 164 @asaplog_post_dec 165 def set_shift(self, mode="DSB", imageshift=None): 166 """ 167 Set shift mode and channel shift of image band. 168 169 mode : shift mode ['DSB'|'SSB'(='2SB')] 170 When mode='DSB', imageshift is assumed to be equal 171 to the shift of signal sideband but in opposite direction. 172 imageshift : a list of number of channel shift in image band of 173 each scantable. valid only mode='SSB' 174 """ 175 if mode.upper().startswith("D"): 176 # DSB mode 177 self.dsbmode = True 178 self.imageShift = [] 179 else: 180 if not imageshift: 181 raise ValueError, "Need to set shift value of image sideband" 182 self.dsbmode = False 183 self.imageShift = imageshift 184 asaplog.push("Image sideband shift is set manually: %s" % str(self.imageShift)) 185 186 @asaplog_post_dec 187 def set_both(self, flag=False): 188 """ 189 Resolve both image and signal sideband when True. 190 """ 191 self.getboth = flag 192 if self.getboth: 193 asaplog.push("Both signal and image sidebands are solved and output as separate tables.") 194 else: 195 asaplog.push("Only signal sideband is solved and output as an table.") 196 197 @asaplog_post_dec 198 def set_limit(self, threshold=0.2): 199 """ 200 Set rejection limit of solution. 201 """ 202 #self._separator._setlimit(abs(threshold)) 203 self.rejlimit = threshold 204 asaplog.push("The threshold of rejection is set to %f" % self.rejlimit) 205 206 207 @asaplog_post_dec 208 def set_solve_other(self, flag=False): 209 """ 210 Calculate spectra by subtracting the solution of the other sideband 211 when True. 212 """ 213 self.solveother = flag 214 if flag: 215 asaplog.push("Expert mode: solution are obtained by subtraction of the other sideband.") 216 217 def set_lo1(self,lo1): 218 """ 219 Set LO1 frequency to calculate frequency of image sideband. 220 221 lo1 : LO1 frequency in float 222 """ 223 lo1val = -1. 224 if isinstance(lo1, dict) and lo1["unit"] == "Hz": 225 lo1val = lo1["value"] 226 else: 227 lo1val = float(lo1) 228 if lo1val <= 0.: 229 asaplog.push("Got negative LO1 frequency. It will be ignored.") 230 asaplog.post("WARN") 231 else: 232 self._separator.set_lo1(lo1val) 233 234 235 def set_lo1root(self, name): 236 """ 237 Set MS name which stores LO1 frequency of signal side band. 238 It is used to calculate frequency of image sideband. 239 240 name : MS name which contains 'ASDM_SPECTRALWINDOW' and 241 'ASDM_RECEIVER' tables. 242 """ 243 self._separator.set_lo1root(name) 244 245 246 @asaplog_post_dec 247 def separate(self, outname="", overwrite=False): 248 """ 249 Invoke sideband separation. 250 251 outname : a name of output scantable 252 overwrite : overwrite existing table 253 """ 254 # List up valid scantables and IFNOs to convolve. 255 #self._separator._separate() 256 self._setup_shift() 257 #self._preprocess_tables() 258 259 ### TEMPORAL ### 260 self._separator._get_asistb_from_scantb(self.tables[0]) 261 ################ 262 263 nshift = len(self.tables) 264 signaltab = self._grid_outtable(self.tables[0]) 265 if self.getboth: 266 imagetab = signaltab.copy() 267 268 rejrow = [] 269 for irow in xrange(signaltab.nrow()): 270 currpol = signaltab.getpol(irow) 271 currbeam = signaltab.getbeam(irow) 272 currdir = signaltab.get_directionval(irow) 273 spec_array, tabidx = self._get_specarray(polid=currpol,\ 274 beamid=currbeam,\ 275 dir=currdir) 276 #if not spec_array: 277 if len(tabidx) == 0: 278 asaplog.post() 279 asaplog.push("skipping row %d" % irow) 280 rejrow.append(irow) 281 continue 282 signal = self._solve_signal(spec_array, tabidx) 283 signaltab.set_spectrum(signal, irow) 284 285 # Solve image side side band 286 if self.getboth: 287 image = self._solve_image(spec_array, tabidx) 288 imagetab.set_spectrum(image, irow) 289 290 # TODO: Need to remove rejrow form scantables here 291 signaltab.flag_row(rejrow) 292 if self.getboth: 293 imagetab.flag_row(rejrow) 169 return False 170 171 self._separator.separate(outname) 172 173 ###################################################################### 174 # @asaplog_post_dec 175 # def separate0(self, outname="", overwrite=False): 176 # """ 177 # Invoke sideband separation. 178 179 # outname : a name of output scantable 180 # overwrite : overwrite existing table 181 # """ 182 # # List up valid scantables and IFNOs to convolve. 183 # #self._separator.separate() 184 # self._setup_shift() 185 # #self._preprocess_tables() 186 187 # ### TEMPORAL ### 188 # self._separator._get_asistb_from_scantb(self.tables[0]) 189 # ################ 190 191 # nshift = len(self.tables) 192 # signaltab = self._grid_outtable(self.tables[0]) 193 # if self.getboth: 194 # imagetab = signaltab.copy() 195 196 # rejrow = [] 197 # for irow in xrange(signaltab.nrow()): 198 # currpol = signaltab.getpol(irow) 199 # currbeam = signaltab.getbeam(irow) 200 # currdir = signaltab.get_directionval(irow) 201 # spec_array, tabidx = self._get_specarray(polid=currpol,\ 202 # beamid=currbeam,\ 203 # dir=currdir) 204 # #if not spec_array: 205 # if len(tabidx) == 0: 206 # asaplog.post() 207 # asaplog.push("skipping row %d" % irow) 208 # rejrow.append(irow) 209 # continue 210 # signal = self._solve_signal(spec_array, tabidx) 211 # signaltab.set_spectrum(signal, irow) 212 213 # # Solve image side side band 214 # if self.getboth: 215 # image = self._solve_image(spec_array, tabidx) 216 # imagetab.set_spectrum(image, irow) 217 218 # # TODO: Need to remove rejrow form scantables here 219 # signaltab.flag_row(rejrow) 220 # if self.getboth: 221 # imagetab.flag_row(rejrow) 294 222 295 if outname == "":296 outname = "sbsepareted.asap"297 signalname = outname + ".signalband"298 if os.path.exists(signalname):299 if not overwrite:300 raise Exception, "Output file '%s' exists." % signalname301 else:302 shutil.rmtree(signalname)303 signaltab.save(signalname)304 305 if self.getboth:306 # Warnings307 asaplog.post()308 asaplog.push("Saving IMAGE sideband.")309 #asaplog.push("Note, frequency information of IMAGE sideband cannot be properly filled so far. (future development)")310 #asaplog.push("Storing frequency setting of SIGNAL sideband in FREQUENCIES table for now.")311 #asaplog.post("WARN")312 313 imagename = outname + ".imageband"314 if os.path.exists(imagename):315 if not overwrite:316 raise Exception, "Output file '%s' exists." % imagename317 else:318 shutil.rmtree(imagename)319 # Update frequency information320 self._separator.set_imgtable(imagetab)321 self._separator.solve_imgfreq()322 imagetab.save(imagename)323 324 325 def _solve_signal(self, data, tabidx=None):326 if not tabidx:327 tabidx = range(len(data))328 329 tempshift = []330 dshift = []331 if self.solveother:332 for idx in tabidx:333 tempshift.append(-self.imageShift[idx])334 dshift.append(self.signalShift[idx] - self.imageShift[idx])335 else:336 for idx in tabidx:337 tempshift.append(-self.signalShift[idx])338 dshift.append(self.imageShift[idx] - self.signalShift[idx])339 340 shiftdata = numpy.zeros(data.shape, numpy.float)341 for i in range(len(data)):342 shiftdata[i] = self._shiftSpectrum(data[i], tempshift[i])343 ifftdata = self._Deconvolution(shiftdata, dshift, self.rejlimit)344 result_image = self._combineResult(ifftdata)345 if not self.solveother:346 return result_image347 result_signal = self._subtractOtherSide(shiftdata, dshift, result_image)348 return result_signal349 350 351 def _solve_image(self, data, tabidx=None):352 if not tabidx:353 tabidx = range(len(data))354 355 tempshift = []356 dshift = []357 if self.solveother:358 for idx in tabidx:359 tempshift.append(-self.signalShift[idx])360 dshift.append(self.imageShift[idx] - self.signalShift[idx])361 else:362 for idx in tabidx:363 tempshift.append(-self.imageShift[idx])364 dshift.append(self.signalShift[idx] - self.imageShift[idx])365 366 shiftdata = numpy.zeros(data.shape, numpy.float)367 for i in range(len(data)):368 shiftdata[i] = self._shiftSpectrum(data[i], tempshift[i])369 ifftdata = self._Deconvolution(shiftdata, dshift, self.rejlimit)370 result_image = self._combineResult(ifftdata)371 if not self.solveother:372 return result_image373 result_signal = self._subtractOtherSide(shiftdata, dshift, result_image)374 return result_signal375 376 @asaplog_post_dec377 def _grid_outtable(self, table):378 # Generate gridded table for output (Just to get rows)379 gridder = asapgrid2(table)380 gridder.setIF(self.baseif)381 382 cellx = str(self.dirtol[0])+"rad"383 celly = str(self.dirtol[1])+"rad"384 dirarr = numpy.array(table.get_directionval()).transpose()385 mapx = dirarr[0].max() - dirarr[0].min()386 mapy = dirarr[1].max() - dirarr[1].min()387 centy = 0.5 * (dirarr[1].max() + dirarr[1].min())388 nx = max(1, numpy.ceil(mapx*numpy.cos(centy)/self.dirtol[0]))389 ny = max(1, numpy.ceil(mapy/self.dirtol[0]))390 391 asaplog.push("Regrid output scantable with cell = [%s, %s]" % \392 (cellx, celly))393 gridder.defineImage(nx=nx, ny=ny, cellx=cellx, celly=celly)394 gridder.setFunc(func='box', convsupport=1)395 gridder.setWeight(weightType='uniform')396 gridder.grid()397 return gridder.getResult()398 399 @asaplog_post_dec400 def _get_specarray(self, polid=None, beamid=None, dir=None):401 ntable = len(self.tables)402 spec_array = numpy.zeros((ntable, self.nchan), numpy.float)403 nspec = 0404 asaplog.push("Start data selection by POL=%d, BEAM=%d, direction=[%f, %f]" % (polid, beamid, dir[0], dir[1]))405 tabidx = []406 for itab in range(ntable):407 tab = self.tables[itab]408 # Select rows by POLNO and BEAMNO409 try:410 tab.set_selection(pols=[polid], beams=[beamid])411 if tab.nrow() > 0: tabidx.append(itab)412 except: # no selection413 asaplog.post()414 asaplog.push("table %d - No spectrum ....skipping the table" % (itab))415 asaplog.post("WARN")416 continue417 418 # Select rows by direction419 spec = numpy.zeros(self.nchan, numpy.float)420 selrow = []421 for irow in range(tab.nrow()):422 currdir = tab.get_directionval(irow)423 if (abs(currdir[0] - dir[0]) > self.dirtol[0]) or \424 (abs(currdir[1] - dir[1]) > self.dirtol[1]):425 continue426 selrow.append(irow)427 if len(selrow) == 0:428 asaplog.post()429 asaplog.push("table %d - No spectrum ....skipping the table" % (itab))430 asaplog.post("WARN")431 continue432 else:433 seltab = tab.copy()434 seltab.set_selection(selector(rows=selrow))435 436 iftab.nrow() > 1:437 asaplog.push("table %d - More than a spectrum selected. averaging rows..." % (itab))438 tab = seltab.average_time(scanav=False, weight="tintsys")439 else:440 tab = seltab441 442 spec_array[nspec] = tab._getspectrum()443 nspec += 1444 445 if nspec != ntable:446 asaplog.post()447 #asaplog.push("Some tables has no spectrum with POL=%d BEAM=%d. averaging rows..." % (polid, beamid))448 asaplog.push("Could not find corresponding rows in some tables.")449 asaplog.push("Number of spectra selected = %d (table: %d)" % (nspec, ntable))450 if nspec < 2:451 asaplog.push("At least 2 spectra are necessary for convolution")452 asaplog.post("ERROR")453 return False, tabidx454 455 return spec_array[0:nspec], tabidx456 457 458 @asaplog_post_dec459 def _setup_shift(self):460 ### define self.tables, self.signalShift, and self.imageShift461 if not self.intables:462 asaplog.post()463 raise RuntimeError, "Input data is not defined."464 #if self.baseif < 0:465 # asaplog.post()466 # raise RuntimeError, "Reference IFNO is not defined."223 # if outname == "": 224 # outname = "sbsepareted.asap" 225 # signalname = outname + ".signalband" 226 # if os.path.exists(signalname): 227 # if not overwrite: 228 # raise Exception, "Output file '%s' exists." % signalname 229 # else: 230 # shutil.rmtree(signalname) 231 # signaltab.save(signalname) 232 233 # if self.getboth: 234 # # Warnings 235 # asaplog.post() 236 # asaplog.push("Saving IMAGE sideband.") 237 # #asaplog.push("Note, frequency information of IMAGE sideband cannot be properly filled so far. (future development)") 238 # #asaplog.push("Storing frequency setting of SIGNAL sideband in FREQUENCIES table for now.") 239 # #asaplog.post("WARN") 240 241 # imagename = outname + ".imageband" 242 # if os.path.exists(imagename): 243 # if not overwrite: 244 # raise Exception, "Output file '%s' exists." % imagename 245 # else: 246 # shutil.rmtree(imagename) 247 # # Update frequency information 248 # self._separator.set_imgtable(imagetab) 249 # self._separator.solve_imgfreq() 250 # imagetab.save(imagename) 251 252 253 # def _solve_signal(self, data, tabidx=None): 254 # if not tabidx: 255 # tabidx = range(len(data)) 256 257 # tempshift = [] 258 # dshift = [] 259 # if self.solveother: 260 # for idx in tabidx: 261 # tempshift.append(-self.imageShift[idx]) 262 # dshift.append(self.signalShift[idx] - self.imageShift[idx]) 263 # else: 264 # for idx in tabidx: 265 # tempshift.append(-self.signalShift[idx]) 266 # dshift.append(self.imageShift[idx] - self.signalShift[idx]) 267 268 # shiftdata = numpy.zeros(data.shape, numpy.float) 269 # for i in range(len(data)): 270 # shiftdata[i] = self._shiftSpectrum(data[i], tempshift[i]) 271 # ifftdata = self._Deconvolution(shiftdata, dshift, self.rejlimit) 272 # result_image = self._combineResult(ifftdata) 273 # if not self.solveother: 274 # return result_image 275 # result_signal = self._subtractOtherSide(shiftdata, dshift, result_image) 276 # return result_signal 277 278 279 # def _solve_image(self, data, tabidx=None): 280 # if not tabidx: 281 # tabidx = range(len(data)) 282 283 # tempshift = [] 284 # dshift = [] 285 # if self.solveother: 286 # for idx in tabidx: 287 # tempshift.append(-self.signalShift[idx]) 288 # dshift.append(self.imageShift[idx] - self.signalShift[idx]) 289 # else: 290 # for idx in tabidx: 291 # tempshift.append(-self.imageShift[idx]) 292 # dshift.append(self.signalShift[idx] - self.imageShift[idx]) 293 294 # shiftdata = numpy.zeros(data.shape, numpy.float) 295 # for i in range(len(data)): 296 # shiftdata[i] = self._shiftSpectrum(data[i], tempshift[i]) 297 # ifftdata = self._Deconvolution(shiftdata, dshift, self.rejlimit) 298 # result_image = self._combineResult(ifftdata) 299 # if not self.solveother: 300 # return result_image 301 # result_signal = self._subtractOtherSide(shiftdata, dshift, result_image) 302 # return result_signal 303 304 # @asaplog_post_dec 305 # def _grid_outtable(self, table): 306 # # Generate gridded table for output (Just to get rows) 307 # gridder = asapgrid2(table) 308 # gridder.setIF(self.baseif) 309 310 # cellx = str(self.dirtol[0])+"rad" 311 # celly = str(self.dirtol[1])+"rad" 312 # dirarr = numpy.array(table.get_directionval()).transpose() 313 # mapx = dirarr[0].max() - dirarr[0].min() 314 # mapy = dirarr[1].max() - dirarr[1].min() 315 # centy = 0.5 * (dirarr[1].max() + dirarr[1].min()) 316 # nx = max(1, numpy.ceil(mapx*numpy.cos(centy)/self.dirtol[0])) 317 # ny = max(1, numpy.ceil(mapy/self.dirtol[0])) 318 319 # asaplog.push("Regrid output scantable with cell = [%s, %s]" % \ 320 # (cellx, celly)) 321 # gridder.defineImage(nx=nx, ny=ny, cellx=cellx, celly=celly) 322 # gridder.setFunc(func='box', convsupport=1) 323 # gridder.setWeight(weightType='uniform') 324 # gridder.grid() 325 # return gridder.getResult() 326 327 # @asaplog_post_dec 328 # def _get_specarray(self, polid=None, beamid=None, dir=None): 329 # ntable = len(self.tables) 330 # spec_array = numpy.zeros((ntable, self.nchan), numpy.float) 331 # nspec = 0 332 # asaplog.push("Start data selection by POL=%d, BEAM=%d, direction=[%f, %f]" % (polid, beamid, dir[0], dir[1])) 333 # tabidx = [] 334 # for itab in range(ntable): 335 # tab = self.tables[itab] 336 # # Select rows by POLNO and BEAMNO 337 # try: 338 # tab.set_selection(pols=[polid], beams=[beamid]) 339 # if tab.nrow() > 0: tabidx.append(itab) 340 # except: # no selection 341 # asaplog.post() 342 # asaplog.push("table %d - No spectrum ....skipping the table" % (itab)) 343 # asaplog.post("WARN") 344 # continue 345 346 # # Select rows by direction 347 # spec = numpy.zeros(self.nchan, numpy.float) 348 # selrow = [] 349 # for irow in range(tab.nrow()): 350 # currdir = tab.get_directionval(irow) 351 # if (abs(currdir[0] - dir[0]) > self.dirtol[0]) or \ 352 # (abs(currdir[1] - dir[1]) > self.dirtol[1]): 353 # continue 354 # selrow.append(irow) 355 # if len(selrow) == 0: 356 # asaplog.post() 357 # asaplog.push("table %d - No spectrum ....skipping the table" % (itab)) 358 # asaplog.post("WARN") 359 # continue 360 # else: 361 # seltab = tab.copy() 362 # seltab.set_selection(selector(rows=selrow)) 363 364 # if seltab.nrow() > 1: 365 # asaplog.push("table %d - More than a spectrum selected. averaging rows..." % (itab)) 366 # tab = seltab.average_time(scanav=False, weight="tintsys") 367 # else: 368 # tab = seltab 369 370 # spec_array[nspec] = tab._getspectrum() 371 # nspec += 1 372 373 # if nspec != ntable: 374 # asaplog.post() 375 # #asaplog.push("Some tables has no spectrum with POL=%d BEAM=%d. averaging rows..." % (polid, beamid)) 376 # asaplog.push("Could not find corresponding rows in some tables.") 377 # asaplog.push("Number of spectra selected = %d (table: %d)" % (nspec, ntable)) 378 # if nspec < 2: 379 # asaplog.push("At least 2 spectra are necessary for convolution") 380 # asaplog.post("ERROR") 381 # return False, tabidx 382 383 # return spec_array[0:nspec], tabidx 384 385 386 # @asaplog_post_dec 387 # def _setup_shift(self): 388 # ### define self.tables, self.signalShift, and self.imageShift 389 # if not self.intables: 390 # asaplog.post() 391 # raise RuntimeError, "Input data is not defined." 392 # #if self.baseif < 0: 393 # # asaplog.post() 394 # # raise RuntimeError, "Reference IFNO is not defined." 467 395 468 byname = False469 #if not self.intables:470 if isinstance(self.intables[0], str):471 # A list of file name is given472 if not os.path.exists(self.intables[0]):473 asaplog.post()474 raise RuntimeError, "Could not find '%s'" % self.intables[0]396 # byname = False 397 # #if not self.intables: 398 # if isinstance(self.intables[0], str): 399 # # A list of file name is given 400 # if not os.path.exists(self.intables[0]): 401 # asaplog.post() 402 # raise RuntimeError, "Could not find '%s'" % self.intables[0] 475 403 476 stab = scantable(self.intables[0],average=False) 477 ntab = len(self.intables) 478 byname = True 479 else: 480 stab = self.intables[0] 481 ntab = len(self.intables) 482 483 if len(stab.getbeamnos()) > 1: 484 asaplog.post() 485 asaplog.push("Mult-beam data is not supported by this module.") 486 asaplog.post("ERROR") 487 return 488 489 valid_ifs = stab.getifnos() 490 if self.baseif < 0: 491 self.baseif = valid_ifs[0] 492 asaplog.post() 493 asaplog.push("IFNO is not selected. Using the first IF in the first scantable. Reference IFNO = %d" % (self.baseif)) 494 495 if not (self.baseif in valid_ifs): 496 asaplog.post() 497 errmsg = "IF%d does not exist in the first scantable" % \ 498 self.baseif 499 raise RuntimeError, errmsg 500 501 asaplog.push("Start selecting tables and IFNOs to solve.") 502 asaplog.push("Checking frequency of the reference IF") 503 unit_org = stab.get_unit() 504 coord = stab._getcoordinfo() 505 frame_org = coord[1] 506 stab.set_unit("Hz") 507 if len(self.freqframe) > 0: 508 stab.set_freqframe(self.freqframe) 509 stab.set_selection(ifs=[self.baseif]) 510 spx = stab._getabcissa() 511 stab.set_selection() 512 basech0 = spx[0] 513 baseinc = spx[1]-spx[0] 514 self.nchan = len(spx) 515 # frequency tolerance 516 if isinstance(self.freqtol, dict) and self.freqtol['unit'] == "Hz": 517 vftol = abs(self.freqtol['value']) 518 else: 519 vftol = abs(baseinc * float(self.freqtol)) 520 self.freqtol = dict(value=vftol, unit="Hz") 521 # tolerance of frequency increment 522 inctol = abs(baseinc/float(self.nchan)) 523 asaplog.push("Reference frequency setup (Table = 0, IFNO = %d): nchan = %d, chan0 = %f Hz, incr = %f Hz" % (self.baseif, self.nchan, basech0, baseinc)) 524 asaplog.push("Allowed frequency tolerance = %f Hz ( %f channels)" % (vftol, vftol/baseinc)) 525 poltype0 = stab.poltype() 526 527 self.tables = [] 528 self.signalShift = [] 529 if self.dsbmode: 530 self.imageShift = [] 531 532 for itab in range(ntab): 533 asaplog.push("Table %d:" % itab) 534 tab_selected = False 535 if itab > 0: 536 if byname: 537 stab = scantable(self.intables[itab],average=False) 538 else: 539 stab = self.intables[itab] 540 unit_org = stab.get_unit() 541 coord = stab._getcoordinfo() 542 frame_org = coord[1] 543 stab.set_unit("Hz") 544 if len(self.freqframe) > 0: 545 stab.set_freqframe(self.freqframe) 546 547 # Check POLTYPE should be equal to the first table. 548 if stab.poltype() != poltype0: 549 asaplog.post() 550 raise Exception, "POLTYPE should be equal to the first table." 551 # Multiple beam data may not handled properly 552 if len(stab.getbeamnos()) > 1: 553 asaplog.post() 554 asaplog.push("table contains multiple beams. It may not be handled properly.") 555 asaplog.push("WARN") 556 557 for ifno in stab.getifnos(): 558 stab.set_selection(ifs=[ifno]) 559 spx = stab._getabcissa() 560 if (len(spx) != self.nchan) or \ 561 (abs(spx[0]-basech0) > vftol) or \ 562 (abs(spx[1]-spx[0]-baseinc) > inctol): 563 continue 564 tab_selected = True 565 seltab = stab.copy() 566 seltab.set_unit(unit_org) 567 seltab.set_freqframe(frame_org) 568 self.tables.append(seltab) 569 self.signalShift.append((spx[0]-basech0)/baseinc) 570 if self.dsbmode: 571 self.imageShift.append(-self.signalShift[-1]) 572 asaplog.push("- IF%d selected: sideband shift = %16.12e channels" % (ifno, self.signalShift[-1])) 573 stab.set_selection() 574 stab.set_unit(unit_org) 575 stab.set_freqframe(frame_org) 576 if not tab_selected: 577 asaplog.post() 578 asaplog.push("No data selected in Table %d" % itab) 579 asaplog.post("WARN") 580 581 asaplog.push("Total number of IFs selected = %d" % len(self.tables)) 582 if len(self.tables) < 2: 583 asaplog.post() 584 raise RuntimeError, "At least 2 IFs are necessary for convolution!" 585 586 if not self.dsbmode and len(self.imageShift) != len(self.signalShift): 587 asaplog.post() 588 errmsg = "User defined channel shift of image sideband has %d elements, while selected IFNOs are %d" % (len(self.imageShift), len(self.signalShift)) 589 errmsg += "\nThe frequency tolerance (freqtol) you set may be too small." 590 raise RuntimeError, errmsg 591 592 self.signalShift = numpy.array(self.signalShift) 593 self.imageShift = numpy.array(self.imageShift) 594 self.nshift = len(self.tables) 595 596 @asaplog_post_dec 597 def _preprocess_tables(self): 598 ### temporary method to preprocess data 599 ### Do time averaging for now. 600 for itab in range(len(self.tables)): 601 self.tables[itab] = self.tables[itab].average_time(scanav=False, weight="tintsys") 404 # stab = scantable(self.intables[0],average=False) 405 # ntab = len(self.intables) 406 # byname = True 407 # else: 408 # stab = self.intables[0] 409 # ntab = len(self.intables) 410 411 # if len(stab.getbeamnos()) > 1: 412 # asaplog.post() 413 # asaplog.push("Mult-beam data is not supported by this module.") 414 # asaplog.post("ERROR") 415 # return 416 417 # valid_ifs = stab.getifnos() 418 # if self.baseif < 0: 419 # self.baseif = valid_ifs[0] 420 # asaplog.post() 421 # asaplog.push("IFNO is not selected. Using the first IF in the first scantable. Reference IFNO = %d" % (self.baseif)) 422 423 # if not (self.baseif in valid_ifs): 424 # asaplog.post() 425 # errmsg = "IF%d does not exist in the first scantable" % \ 426 # self.baseif 427 # raise RuntimeError, errmsg 428 429 # asaplog.push("Start selecting tables and IFNOs to solve.") 430 # asaplog.push("Checking frequency of the reference IF") 431 # unit_org = stab.get_unit() 432 # coord = stab._getcoordinfo() 433 # frame_org = coord[1] 434 # stab.set_unit("Hz") 435 # if len(self.freqframe) > 0: 436 # stab.set_freqframe(self.freqframe) 437 # stab.set_selection(ifs=[self.baseif]) 438 # spx = stab._getabcissa() 439 # stab.set_selection() 440 # basech0 = spx[0] 441 # baseinc = spx[1]-spx[0] 442 # self.nchan = len(spx) 443 # # frequency tolerance 444 # if isinstance(self.freqtol, dict) and self.freqtol['unit'] == "Hz": 445 # vftol = abs(self.freqtol['value']) 446 # else: 447 # vftol = abs(baseinc * float(self.freqtol)) 448 # self.freqtol = dict(value=vftol, unit="Hz") 449 # # tolerance of frequency increment 450 # inctol = abs(baseinc/float(self.nchan)) 451 # asaplog.push("Reference frequency setup (Table = 0, IFNO = %d): nchan = %d, chan0 = %f Hz, incr = %f Hz" % (self.baseif, self.nchan, basech0, baseinc)) 452 # asaplog.push("Allowed frequency tolerance = %f Hz ( %f channels)" % (vftol, vftol/baseinc)) 453 # poltype0 = stab.poltype() 454 455 # self.tables = [] 456 # self.signalShift = [] 457 # if self.dsbmode: 458 # self.imageShift = [] 459 460 # for itab in range(ntab): 461 # asaplog.push("Table %d:" % itab) 462 # tab_selected = False 463 # if itab > 0: 464 # if byname: 465 # stab = scantable(self.intables[itab],average=False) 466 # else: 467 # stab = self.intables[itab] 468 # unit_org = stab.get_unit() 469 # coord = stab._getcoordinfo() 470 # frame_org = coord[1] 471 # stab.set_unit("Hz") 472 # if len(self.freqframe) > 0: 473 # stab.set_freqframe(self.freqframe) 474 475 # # Check POLTYPE should be equal to the first table. 476 # if stab.poltype() != poltype0: 477 # asaplog.post() 478 # raise Exception, "POLTYPE should be equal to the first table." 479 # # Multiple beam data may not handled properly 480 # if len(stab.getbeamnos()) > 1: 481 # asaplog.post() 482 # asaplog.push("table contains multiple beams. It may not be handled properly.") 483 # asaplog.push("WARN") 484 485 # for ifno in stab.getifnos(): 486 # stab.set_selection(ifs=[ifno]) 487 # spx = stab._getabcissa() 488 # if (len(spx) != self.nchan) or \ 489 # (abs(spx[0]-basech0) > vftol) or \ 490 # (abs(spx[1]-spx[0]-baseinc) > inctol): 491 # continue 492 # tab_selected = True 493 # seltab = stab.copy() 494 # seltab.set_unit(unit_org) 495 # seltab.set_freqframe(frame_org) 496 # self.tables.append(seltab) 497 # self.signalShift.append((spx[0]-basech0)/baseinc) 498 # if self.dsbmode: 499 # self.imageShift.append(-self.signalShift[-1]) 500 # asaplog.push("- IF%d selected: sideband shift = %16.12e channels" % (ifno, self.signalShift[-1])) 501 # stab.set_selection() 502 # stab.set_unit(unit_org) 503 # stab.set_freqframe(frame_org) 504 # if not tab_selected: 505 # asaplog.post() 506 # asaplog.push("No data selected in Table %d" % itab) 507 # asaplog.post("WARN") 508 509 # asaplog.push("Total number of IFs selected = %d" % len(self.tables)) 510 # if len(self.tables) < 2: 511 # asaplog.post() 512 # raise RuntimeError, "At least 2 IFs are necessary for convolution!" 513 514 # if not self.dsbmode and len(self.imageShift) != len(self.signalShift): 515 # asaplog.post() 516 # errmsg = "User defined channel shift of image sideband has %d elements, while selected IFNOs are %d" % (len(self.imageShift), len(self.signalShift)) 517 # errmsg += "\nThe frequency tolerance (freqtol) you set may be too small." 518 # raise RuntimeError, errmsg 519 520 # self.signalShift = numpy.array(self.signalShift) 521 # self.imageShift = numpy.array(self.imageShift) 522 # self.nshift = len(self.tables) 602 523 603 524 # @asaplog_post_dec 604 # def _setup_image_freq(self, table): 605 # # User defined coordinate 606 # # Get from associated MS 607 # # Get from ASDM 608 # lo1 = -1. 609 # if self.lo1 > 0.: 610 # asaplog.push("Using user defined LO1 frequency %e16.12 [Hz]" % self.lo1) 611 # lo1 = self.lo1 612 # else: 613 # print "NOT IMPLEMENTED YET!!!" 614 615 # def save(self, outfile, outform="ASAP", overwrite=False): 616 # if not overwrite and os.path.exists(outfile): 617 # raise RuntimeError, "Output file '%s' already exists" % outfile 618 # 619 # #self._separator._save(outfile, outform) 620 621 # def done(self): 622 # self.close() 623 624 # def close(self): 625 # pass 626 # #del self._separator 627 628 629 630 ######################################################################## 631 def _Deconvolution(self, data_array, shift, threshold=0.00000001): 632 FObs = [] 633 Reject = 0 634 nshift, nchan = data_array.shape 635 nspec = nshift*(nshift-1)/2 636 ifftObs = numpy.zeros((nspec, nchan), numpy.float) 637 for i in range(nshift): 638 F = FFT.fft(data_array[i]) 639 FObs.append(F) 640 z = 0 641 for i in range(nshift): 642 for j in range(i+1, nshift): 643 Fobs = (FObs[i]+FObs[j])/2.0 644 dX = (shift[j]-shift[i])*2.0*math.pi/float(self.nchan) 645 #print 'dX,i,j=',dX,i,j 646 for k in range(1,self.nchan): 647 if math.fabs(math.sin(dX*k)) > threshold: 648 Fobs[k] += ((FObs[i][k]-FObs[j][k])/2.0/(1.0-math.cos(dX*k))*math.sin(dX*k))*1.0j 649 else: Reject += 1 650 ifftObs[z] = FFT.ifft(Fobs) 651 z += 1 652 print 'Threshold=%s Reject=%d' % (threshold, Reject) 653 return ifftObs 654 655 def _combineResult(self, ifftObs): 656 nspec = len(ifftObs) 657 sum = ifftObs[0] 658 for i in range(1,nspec): 659 sum += ifftObs[i] 660 return(sum/float(nspec)) 661 662 def _subtractOtherSide(self, data_array, shift, Data): 663 sum = numpy.zeros(len(Data), numpy.float) 664 numSP = len(data_array) 665 for i in range(numSP): 666 SPsub = data_array[i] - Data 667 sum += self._shiftSpectrum(SPsub, -shift[i]) 668 return(sum/float(numSP)) 669 670 def _shiftSpectrum(self, data, Shift): 671 Out = numpy.zeros(self.nchan, numpy.float) 672 w2 = Shift % 1 673 w1 = 1.0 - w2 674 for i in range(self.nchan): 675 c1 = int((Shift + i) % self.nchan) 676 c2 = (c1 + 1) % self.nchan 677 Out[c1] += data[i] * w1 678 Out[c2] += data[i] * w2 679 return Out.copy() 525 # def _preprocess_tables(self): 526 # ### temporary method to preprocess data 527 # ### Do time averaging for now. 528 # for itab in range(len(self.tables)): 529 # self.tables[itab] = self.tables[itab].average_time(scanav=False, weight="tintsys") 530 531 532 # ######################################################################## 533 # def _Deconvolution(self, data_array, shift, threshold=0.00000001): 534 # FObs = [] 535 # Reject = 0 536 # nshift, nchan = data_array.shape 537 # nspec = nshift*(nshift-1)/2 538 # ifftObs = numpy.zeros((nspec, nchan), numpy.float) 539 # for i in range(nshift): 540 # F = FFT.fft(data_array[i]) 541 # FObs.append(F) 542 # z = 0 543 # for i in range(nshift): 544 # for j in range(i+1, nshift): 545 # Fobs = (FObs[i]+FObs[j])/2.0 546 # dX = (shift[j]-shift[i])*2.0*math.pi/float(self.nchan) 547 # #print 'dX,i,j=',dX,i,j 548 # for k in range(1,self.nchan): 549 # if math.fabs(math.sin(dX*k)) > threshold: 550 # Fobs[k] += ((FObs[i][k]-FObs[j][k])/2.0/(1.0-math.cos(dX*k))*math.sin(dX*k))*1.0j 551 # else: Reject += 1 552 # ifftObs[z] = FFT.ifft(Fobs) 553 # z += 1 554 # print 'Threshold=%s Reject=%d' % (threshold, Reject) 555 # return ifftObs 556 557 # def _combineResult(self, ifftObs): 558 # nspec = len(ifftObs) 559 # sum = ifftObs[0] 560 # for i in range(1,nspec): 561 # sum += ifftObs[i] 562 # return(sum/float(nspec)) 563 564 # def _subtractOtherSide(self, data_array, shift, Data): 565 # sum = numpy.zeros(len(Data), numpy.float) 566 # numSP = len(data_array) 567 # for i in range(numSP): 568 # SPsub = data_array[i] - Data 569 # sum += self._shiftSpectrum(SPsub, -shift[i]) 570 # return(sum/float(numSP)) 571 572 # def _shiftSpectrum(self, data, Shift): 573 # Out = numpy.zeros(self.nchan, numpy.float) 574 # w2 = Shift % 1 575 # w1 = 1.0 - w2 576 # for i in range(self.nchan): 577 # c1 = int((Shift + i) % self.nchan) 578 # c2 = (c1 + 1) % self.nchan 579 # Out[c1] += data[i] * w1 580 # Out[c2] += data[i] * w2 581 # return Out.copy() -
trunk/src/STSideBandSep.cpp
r2726 r2794 26 26 27 27 // asap 28 #include "STGrid.h" 29 #include "STMath.h" 30 #include "MathUtils.h" 28 31 #include "STSideBandSep.h" 29 32 … … 32 35 using namespace asap ; 33 36 34 #ifndef KS_DEBUG35 #define KS_DEBUG36 #endif37 // #ifndef KS_DEBUG 38 // #define KS_DEBUG 39 // #endif 37 40 38 41 namespace asap { … … 78 81 79 82 init(); 83 tp_ = intabList_[0]->table().tableType(); 80 84 81 85 os << ntable_ << " tables are set." << LogIO::POST; … … 92 96 { 93 97 // frequency setup 94 sigIfno_= 0;98 sigIfno_= -1; 95 99 ftol_ = -1; 96 100 solFrame_ = MFrequency::N_Types; … … 109 113 // Default LO frame is TOPO 110 114 loFrame_ = MFrequency::TOPO; 115 // scantable storage 116 tp_ = Table::Memory; 111 117 }; 112 118 … … 130 136 131 137 // IFNO 132 sigIfno_ = ifno;138 sigIfno_ = (int) ifno; 133 139 134 140 // Frequency tolerance … … 219 225 if (i != imgShift_.size()-1) os << " , "; 220 226 } 221 os << " ) [channel ]" << LogIO::POST;227 os << " ) [channels]" << LogIO::POST; 222 228 } 223 229 }; … … 227 233 LogIO os(LogOrigin("STSideBandSep","setThreshold()", WHERE)); 228 234 if (limit < 0) 229 throw( AipsError("Rejection limit should be positive number.") );235 throw( AipsError("Rejection limit should be a positive number.") ); 230 236 231 237 rejlimit_ = limit; 232 os << "Rejection limit is set to " << rejlimit_; 233 }; 234 235 // Temporal function to set scantable of image sideband 236 void STSideBandSep::setImageTable(const ScantableWrapper &s) 237 { 238 #ifdef KS_DEBUG 239 cout << "STSideBandSep::setImageTable" << endl; 240 cout << "got image table nrow = " << s.nrow() << endl; 241 #endif 242 imgTab_p = s.getCP(); 243 AlwaysAssert(!imgTab_p.null(),AipsError); 244 245 }; 246 247 // 248 void STSideBandSep::setLO1(const double lo1, const string frame, 238 os << "Rejection limit is set to " << rejlimit_ << LogIO::POST; 239 }; 240 241 void STSideBandSep::separate(string outname) 242 { 243 #ifdef KS_DEBUG 244 cout << "STSideBandSep::separate" << endl; 245 #endif 246 LogIO os(LogOrigin("STSideBandSep","separate()", WHERE)); 247 if (outname.empty()) 248 outname = "sbseparated.asap"; 249 250 // Set up a goup of IFNOs in the list of scantables within 251 // the frequency tolerance and make them a list. 252 nshift_ = setupShift(); 253 if (nshift_ < 2) 254 throw( AipsError("At least 2 IFs are necessary for convolution.") ); 255 // Grid scantable and generate output tables 256 ScantableWrapper gridst = gridTable(); 257 sigTab_p = gridst.getCP(); 258 if (doboth_) 259 imgTab_p = gridst.copy().getCP(); 260 vector<unsigned int> remRowIds; 261 remRowIds.resize(0); 262 Matrix<float> specMat(nchan_, nshift_); 263 vector<float> sigSpec(nchan_), imgSpec(nchan_); 264 vector<uInt> tabIdvec; 265 266 //Generate FFTServer 267 fftsf.resize(IPosition(1, nchan_), FFTEnums::REALTOCOMPLEX); 268 fftsi.resize(IPosition(1, nchan_), FFTEnums::COMPLEXTOREAL); 269 270 /// Loop over sigTab_p and separate sideband 271 for (int irow = 0; irow < sigTab_p->nrow(); irow++){ 272 tabIdvec.resize(0); 273 const int polId = sigTab_p->getPol(irow); 274 const int beamId = sigTab_p->getBeam(irow); 275 const vector<double> dir = sigTab_p->getDirectionVector(irow); 276 // Get a set of spectra to solve 277 if (!getSpectraToSolve(polId, beamId, dir[0], dir[1], 278 specMat, tabIdvec)){ 279 remRowIds.push_back(irow); 280 #ifdef KS_DEBUG 281 cout << "no matching row found. skipping row = " << irow << endl; 282 #endif 283 continue; 284 } 285 // Solve signal sideband 286 sigSpec = solve(specMat, tabIdvec, true); 287 sigTab_p->setSpectrum(sigSpec, irow); 288 289 // Solve image sideband 290 if (doboth_) { 291 imgSpec = solve(specMat, tabIdvec, false); 292 imgTab_p->setSpectrum(imgSpec, irow); 293 } 294 } // end of row loop 295 296 // Remove or flag rows without relevant data from gridded tables 297 if (remRowIds.size() > 0) { 298 const size_t nrem = remRowIds.size(); 299 if ( sigTab_p->table().canRemoveRow() ) { 300 sigTab_p->table().removeRow(remRowIds); 301 os << "Removing " << nrem << " rows from the signal band table" 302 << LogIO::POST; 303 } else { 304 sigTab_p->flagRow(remRowIds, false); 305 os << "Cannot remove rows from the signal band table. Flagging " 306 << nrem << " rows" << LogIO::POST; 307 } 308 309 if (doboth_) { 310 if ( imgTab_p->table().canRemoveRow() ) { 311 imgTab_p->table().removeRow(remRowIds); 312 os << "Removing " << nrem << " rows from the image band table" 313 << LogIO::POST; 314 } else { 315 imgTab_p->flagRow(remRowIds, false); 316 os << "Cannot remove rows from the image band table. Flagging " 317 << nrem << " rows" << LogIO::POST; 318 } 319 } 320 } 321 322 // Finally, save tables on disk 323 if (outname.size() ==0) 324 outname = "sbseparated.asap"; 325 const string sigName = outname + ".signalband"; 326 os << "Saving SIGNAL sideband table: " << sigName << LogIO::POST; 327 sigTab_p->makePersistent(sigName); 328 if (doboth_) { 329 solveImageFrequency(); 330 const string imgName = outname + ".imageband"; 331 os << "Saving IMAGE sideband table: " << sigName << LogIO::POST; 332 imgTab_p->makePersistent(imgName); 333 } 334 335 }; 336 337 unsigned int STSideBandSep::setupShift() 338 { 339 LogIO os(LogOrigin("STSideBandSep","setupShift()", WHERE)); 340 if (infileList_.size() == 0 && intabList_.size() == 0) 341 throw( AipsError("No scantable has been set. Set a list of scantables first.") ); 342 343 const bool byname = (intabList_.size() == 0); 344 // Make sure sigIfno_ exists in the first table. 345 CountedPtr<Scantable> stab; 346 vector<string> coordsav; 347 vector<string> coordtmp(3); 348 os << "Checking IFNO in the first table." << LogIO::POST; 349 if (byname) { 350 if (!checkFile(infileList_[0], "d")) 351 os << LogIO::SEVERE << "Could not find scantable '" << infileList_[0] 352 << "'" << LogIO::EXCEPTION; 353 stab = CountedPtr<Scantable>(new Scantable(infileList_[0])); 354 } else { 355 stab = intabList_[0]; 356 } 357 if (sigIfno_ < 0) { 358 sigIfno_ = (int) stab->getIF(0); 359 os << "IFNO to process has not been set. Using the first IF = " 360 << sigIfno_ << LogIO::POST; 361 } 362 363 unsigned int basench; 364 double basech0, baseinc, ftolval, inctolval; 365 coordsav = stab->getCoordInfo(); 366 const string stfframe = coordsav[1]; 367 coordtmp[0] = "Hz"; 368 coordtmp[1] = ( (solFrame_ == MFrequency::N_Types) ? 369 stfframe : 370 MFrequency::showType(solFrame_) ); 371 coordtmp[2] = coordsav[2]; 372 stab->setCoordInfo(coordtmp); 373 if (!getFreqInfo(stab, (unsigned int) sigIfno_, basech0, baseinc, basench)) { 374 os << LogIO::SEVERE << "No data with IFNO=" << sigIfno_ 375 << " found in the first table." << LogIO::EXCEPTION; 376 } 377 else { 378 os << "Found IFNO = " << sigIfno_ 379 << " in the first table." << LogIO::POST; 380 } 381 if (ftol_.getUnit().empty()) { 382 // tolerance in unit of channels 383 ftolval = ftol_.getValue() * baseinc; 384 } 385 else { 386 ftolval = ftol_.getValue("Hz"); 387 } 388 inctolval = abs(baseinc/(double) basench); 389 const string poltype0 = stab->getPolType(); 390 391 // Initialize shift values 392 initshift(); 393 394 const bool setImg = ( doboth_ && (imgShift_.size() == 0) ); 395 // Select IFs 396 for (unsigned int itab = 0; itab < ntable_; itab++ ){ 397 os << "Table " << itab << LogIO::POST; 398 if (itab > 0) { 399 if (byname) { 400 if (!checkFile(infileList_[itab], "d")) 401 os << LogIO::SEVERE << "Could not find scantable '" 402 << infileList_[itab] << "'" << LogIO::EXCEPTION; 403 stab = CountedPtr<Scantable>(new Scantable(infileList_[itab])); 404 } else { 405 stab = intabList_[itab]; 406 } 407 //POLTYPE should be the same. 408 if (stab->getPolType() != poltype0 ) { 409 os << LogIO::WARN << "POLTYPE differs from the first table." 410 << " Skipping the table" << LogIO::POST; 411 continue; 412 } 413 // Multi beam data may not handled properly 414 if (stab->nbeam() > 1) 415 os << LogIO::WARN << "Table contains multiple beams. " 416 << "It may not be handled properly." << LogIO::POST; 417 418 coordsav = stab->getCoordInfo(); 419 coordtmp[2] = coordsav[2]; 420 stab->setCoordInfo(coordtmp); 421 } 422 bool selected = false; 423 vector<uint> ifnos = stab->getIFNos(); 424 vector<uint>::iterator iter; 425 const STSelector& basesel = stab->getSelection(); 426 for (iter = ifnos.begin(); iter != ifnos.end(); iter++){ 427 unsigned int nch; 428 double freq0, incr; 429 if ( getFreqInfo(stab, *iter, freq0, incr, nch) ){ 430 if ( (nch == basench) && (abs(freq0-basech0) < ftolval) 431 && (abs(incr-baseinc) < inctolval) ){ 432 //Found 433 STSelector sel(basesel); 434 sel.setIFs(vector<int>(1,(int) *iter)); 435 stab->setSelection(sel); 436 CountedPtr<Scantable> seltab = ( new Scantable(*stab, false) ); 437 stab->setSelection(basesel); 438 seltab->setCoordInfo(coordsav); 439 const double chShift = (freq0 - basech0) / baseinc; 440 tableList_.push_back(seltab); 441 sigShift_.push_back(chShift); 442 if (setImg) 443 imgShift_.push_back(-chShift); 444 445 selected = true; 446 os << "- IF" << *iter << " selected: sideband shift = " 447 << chShift << " channels" << LogIO::POST; 448 } 449 } 450 } // ifno loop 451 stab->setCoordInfo(coordsav); 452 if (!selected) 453 os << LogIO::WARN << "No data selected in Table " 454 << itab << LogIO::POST; 455 } // table loop 456 nchan_ = basench; 457 458 os << "Total number of IFs selected = " << tableList_.size() 459 << LogIO::POST; 460 if ( setImg && (sigShift_.size() != imgShift_.size()) ){ 461 os << LogIO::SEVERE 462 << "User defined channel shift of image sideband has " 463 << imgShift_.size() 464 << " elements, while selected IFNOs are " << sigShift_.size() 465 << "\nThe frequency tolerance (freqtol) may be too small." 466 << LogIO::EXCEPTION; 467 } 468 469 return tableList_.size(); 470 }; 471 472 bool STSideBandSep::getFreqInfo(const CountedPtr<Scantable> &stab, 473 const unsigned int &ifno, 474 double &freq0, double &incr, 475 unsigned int &nchan) 476 { 477 vector<uint> ifnos = stab->getIFNos(); 478 bool found = false; 479 vector<uint>::iterator iter; 480 for (iter = ifnos.begin(); iter != ifnos.end(); iter++){ 481 if (*iter == ifno) { 482 found = true; 483 break; 484 } 485 } 486 if (!found) 487 return false; 488 489 const STSelector& basesel = stab->getSelection(); 490 STSelector sel(basesel); 491 sel.setIFs(vector<int>(1,(int) ifno)); 492 stab->setSelection(sel); 493 vector<double> freqs; 494 freqs = stab->getAbcissa(0); 495 freq0 = freqs[0]; 496 incr = freqs[1] - freqs[0]; 497 nchan = freqs.size(); 498 stab->setSelection(basesel); 499 return true; 500 }; 501 502 ScantableWrapper STSideBandSep::gridTable() 503 { 504 LogIO os(LogOrigin("STSideBandSep","gridTable()", WHERE)); 505 if (tableList_.size() == 0) 506 throw( AipsError("Internal error. No scantable has been set to grid.") ); 507 Double xmin, xmax, ymin, ymax; 508 mapExtent(tableList_, xmin, xmax, ymin, ymax); 509 const Double centx = 0.5 * (xmin + xmax); 510 const Double centy = 0.5 * (ymin + ymax); 511 const int nx = max(1, (int) ceil( (xmax - xmin) * cos(centy) /xtol_ ) ); 512 const int ny = max(1, (int) ceil( (ymax - ymin) / ytol_ ) ); 513 514 string scellx, scelly; 515 { 516 ostringstream oss; 517 oss << xtol_ << "rad" ; 518 scellx = oss.str(); 519 } 520 { 521 ostringstream oss; 522 oss << ytol_ << "rad" ; 523 scelly = oss.str(); 524 } 525 526 ScantableWrapper stab0; 527 if (intabList_.size() > 0) 528 stab0 = ScantableWrapper(intabList_[0]); 529 else 530 stab0 = ScantableWrapper(infileList_[0]); 531 532 string scenter; 533 { 534 ostringstream oss; 535 oss << stab0.getCP()->getDirectionRefString() << " " 536 << centx << "rad" << " " << centy << "rad"; 537 scenter = oss.str(); 538 } 539 540 STGrid2 gridder = STGrid2(stab0); 541 gridder.setIF(sigIfno_); 542 gridder.defineImage(nx, ny, scellx, scelly, scenter); 543 gridder.setFunc("box", 1); 544 gridder.setWeight("uniform"); 545 #ifdef KS_DEBUG 546 cout << "Grid parameter summary: " << endl; 547 cout << "- IF = " << sigIfno_ << endl; 548 cout << "- center = " << scenter << ")\n" 549 << "- npix = (" << nx << ", " << ny << ")\n" 550 << "- cell = (" << scellx << ", " << scelly << endl; 551 #endif 552 gridder.grid(); 553 const int itp = (tp_ == Table::Memory ? 0 : 1); 554 ScantableWrapper gtab = gridder.getResultAsScantable(itp); 555 return gtab; 556 }; 557 558 void STSideBandSep::mapExtent(vector< CountedPtr<Scantable> > &tablist, 559 Double &xmin, Double &xmax, 560 Double &ymin, Double &ymax) 561 { 562 ROArrayColumn<Double> dirCol_; 563 dirCol_.attach( tablist[0]->table(), "DIRECTION" ); 564 Matrix<Double> direction = dirCol_.getColumn(); 565 Vector<Double> ra( direction.row(0) ); 566 mathutil::rotateRA(ra); 567 minMax( xmin, xmax, ra ); 568 minMax( ymin, ymax, direction.row(1) ); 569 Double amin, amax, bmin, bmax; 570 const uInt ntab = tablist.size(); 571 for ( uInt i = 1 ; i < ntab ; i++ ) { 572 dirCol_.attach( tablist[i]->table(), "DIRECTION" ); 573 direction.assign( dirCol_.getColumn() ); 574 ra.assign( direction.row(0) ); 575 mathutil::rotateRA(ra); 576 minMax( amin, amax, ra ); 577 minMax( bmin, bmax, direction.row(1) ); 578 xmin = min(xmin, amin); 579 xmax = max(xmax, amax); 580 ymin = min(ymin, bmin); 581 ymax = max(ymax, bmax); 582 } 583 }; 584 585 bool STSideBandSep::getSpectraToSolve(const int polId, const int beamId, 586 const double dirX, const double dirY, 587 Matrix<float> &specMat, vector<uInt> &tabIdvec) 588 { 589 LogIO os(LogOrigin("STSideBandSep","getSpectraToSolve()", WHERE)); 590 591 tabIdvec.resize(0); 592 specMat.resize(nchan_, nshift_); 593 Vector<float> spec; 594 uInt nspec = 0; 595 STMath stm(false); // insitu has no effect for average. 596 for (uInt itab = 0 ; itab < nshift_ ; itab++) { 597 CountedPtr<Scantable> currtab_p = tableList_[itab]; 598 // Selection by POLNO and BEAMNO 599 const STSelector& basesel = currtab_p->getSelection(); 600 STSelector sel(basesel); 601 sel.setPolarizations(vector<int>(1, polId)); 602 sel.setBeams(vector<int>(1, beamId)); 603 try { 604 currtab_p->setSelection(sel); 605 } catch (...) { 606 #ifdef KS_DEBUG 607 cout << "Table " << itab << " - No spectrum found. skipping the table." 608 << endl; 609 #endif 610 continue; 611 } 612 // Selection by direction; 613 vector<int> selrow(0); 614 vector<double> currDir(2, 0.); 615 const int nselrow = currtab_p->nrow(); 616 for (int irow = 0 ; irow < nselrow ; irow++) { 617 currDir = currtab_p->getDirectionVector(irow); 618 if ( (abs(currDir[0]-dirX) > xtol_) || 619 (abs(currDir[1]-dirY) > ytol_) ) 620 continue; 621 // within direction tolerance 622 selrow.push_back(irow); 623 } // end of row loop 624 625 if (selrow.size() < 1) { 626 currtab_p->setSelection(basesel); 627 628 #ifdef KS_DEBUG 629 cout << "Table " << itab << " - No spectrum found. skipping the table." 630 << endl; 631 #endif 632 633 continue; 634 } 635 636 // At least a spectrum is selected in this table 637 CountedPtr<Scantable> seltab_p = ( new Scantable(*currtab_p, false) ); 638 currtab_p->setSelection(basesel); 639 STSelector sel2(seltab_p->getSelection()); 640 sel2.setRows(selrow); 641 seltab_p->setSelection(sel2); 642 CountedPtr<Scantable> avetab_p; 643 vector<bool> mask; 644 if (seltab_p->nrow() > 1) { 645 avetab_p = stm.average(vector< CountedPtr<Scantable> >(1, seltab_p), 646 vector<bool>(), "TINTSYS", "NONE"); 647 #ifdef KS_DEBUG 648 cout << "Table " << itab << " - more than a spectrum is selected. averaging rows..." 649 << endl; 650 #endif 651 if (avetab_p->nrow() > 1) 652 throw( AipsError("Averaged table has more than a row. Somethigs went wrong.") ); 653 } else { 654 avetab_p = seltab_p; 655 } 656 spec.reference(specMat.column(nspec)); 657 spec = avetab_p->getSpectrum(0); 658 tabIdvec.push_back((uInt) itab); 659 nspec++; 660 } // end of table loop 661 if (nspec != nshift_){ 662 //shiftSpecmat.resize(nchan_, nspec, true); 663 #ifdef KS_DEBUG 664 cout << "Could not find corresponding rows in some tables." 665 << endl; 666 cout << "Number of spectra selected = " << nspec << endl; 667 #endif 668 } 669 if (nspec < 2) { 670 #ifdef KS_DEBUG 671 cout << "At least 2 spectra are necessary for convolution" 672 << endl; 673 #endif 674 return false; 675 } 676 return true; 677 }; 678 679 vector<float> STSideBandSep::solve(const Matrix<float> &specmat, 680 const vector<uInt> &tabIdvec, 681 const bool signal) 682 { 683 LogIO os(LogOrigin("STSideBandSep","solve()", WHERE)); 684 if (tabIdvec.size() == 0) 685 throw(AipsError("Internal error. Table index is not defined.")); 686 if (specmat.ncolumn() != tabIdvec.size()) 687 throw(AipsError("Internal error. The row number of input matrix is not conformant.")); 688 if (specmat.nrow() != nchan_) 689 throw(AipsError("Internal error. The channel size of input matrix is not conformant.")); 690 691 692 #ifdef KS_DEBUG 693 cout << "Solving " << (signal ? "SIGNAL" : "IMAGE") << " sideband." 694 << endl; 695 #endif 696 697 const size_t nspec = tabIdvec.size(); 698 vector<double> *thisShift, *otherShift; 699 if (signal == otherside_) { 700 // (solve signal && solveother = T) OR (solve image && solveother = F) 701 thisShift = &imgShift_; 702 otherShift = &sigShift_; 703 #ifdef KS_DEBUG 704 cout << "Image sideband will be deconvolved." << endl; 705 #endif 706 } else { 707 // (solve signal && solveother = F) OR (solve image && solveother = T) 708 thisShift = &sigShift_; 709 otherShift = &imgShift_; 710 #ifdef KS_DEBUG 711 cout << "Signal sideband will be deconvolved." << endl; 712 #endif 713 } 714 715 vector<double> spshift(nspec); 716 Matrix<float> shiftSpecmat(nchan_, nspec, 0.); 717 double tempshift; 718 Vector<float> shiftspvec; 719 for (uInt i = 0 ; i < nspec; i++) { 720 spshift[i] = otherShift->at(i) - thisShift->at(i); 721 tempshift = - thisShift->at(i); 722 shiftspvec.reference(shiftSpecmat.column(i)); 723 shiftSpectrum(specmat.column(i), tempshift, shiftspvec); 724 } 725 726 Matrix<float> convmat(nchan_, nspec*(nspec-1)/2, 0.); 727 vector<float> thisvec(nchan_, 0.); 728 729 float minval, maxval; 730 minMax(minval, maxval, shiftSpecmat); 731 #ifdef KS_DEBUG 732 cout << "Max/Min of input Matrix = (max: " << maxval << ", min: " << minval << ")" << endl; 733 #endif 734 735 #ifdef KS_DEBUG 736 cout << "starting deconvolution" << endl; 737 #endif 738 deconvolve(shiftSpecmat, spshift, rejlimit_, convmat); 739 #ifdef KS_DEBUG 740 cout << "finished deconvolution" << endl; 741 #endif 742 743 minMax(minval, maxval, convmat); 744 #ifdef KS_DEBUG 745 cout << "Max/Min of output Matrix = (max: " << maxval << ", min: " << minval << ")" << endl; 746 #endif 747 748 aggregateMat(convmat, thisvec); 749 750 if (!otherside_) return thisvec; 751 752 // subtract from the other side band. 753 vector<float> othervec(nchan_, 0.); 754 subtractFromOther(shiftSpecmat, thisvec, spshift, othervec); 755 return othervec; 756 }; 757 758 759 void STSideBandSep::shiftSpectrum(const Vector<float> &invec, 760 double shift, 761 Vector<float> &outvec) 762 { 763 LogIO os(LogOrigin("STSideBandSep","shiftSpectrum()", WHERE)); 764 if (invec.size() != nchan_) 765 throw(AipsError("Internal error. The length of input vector differs from nchan_")); 766 if (outvec.size() != nchan_) 767 throw(AipsError("Internal error. The length of output vector differs from nchan_")); 768 769 #ifdef KS_DEBUG 770 cout << "Start shifting spectrum for " << shift << "channels" << endl; 771 #endif 772 773 // tweak shift to be in 0 ~ nchan_-1 774 if ( fabs(shift) > nchan_ ) shift = fmod(shift, nchan_); 775 if (shift < 0.) shift += nchan_; 776 double rweight = fmod(shift, 1.); 777 if (rweight < 0.) rweight += 1.; 778 double lweight = 1. - rweight; 779 uInt lchan, rchan; 780 781 outvec = 0; 782 for (uInt i = 0 ; i < nchan_ ; i++) { 783 lchan = uInt( floor( fmod( (i + shift), nchan_ ) ) ); 784 if (lchan < 0.) lchan += nchan_; 785 rchan = ( (lchan + 1) % nchan_ ); 786 outvec(lchan) += invec(i) * lweight; 787 outvec(rchan) += invec(i) * rweight; 788 } 789 }; 790 791 void STSideBandSep::deconvolve(Matrix<float> &specmat, 792 const vector<double> shiftvec, 793 const double threshold, 794 Matrix<float> &outmat) 795 { 796 LogIO os(LogOrigin("STSideBandSep","deconvolve()", WHERE)); 797 if (specmat.nrow() != nchan_) 798 throw(AipsError("Internal error. The length of input matrix differs from nchan_")); 799 if (specmat.ncolumn() != shiftvec.size()) 800 throw(AipsError("Internal error. The number of input shifts and spectrum differs.")); 801 802 float minval, maxval; 803 #ifdef KS_DEBUG 804 minMax(minval, maxval, specmat); 805 cout << "Max/Min of input Matrix = (max: " << maxval << ", min: " << minval << ")" << endl; 806 #endif 807 808 uInt ninsp = shiftvec.size(); 809 outmat.resize(nchan_, ninsp*(ninsp-1)/2, 0.); 810 Matrix<Complex> fftspmat(nchan_/2+1, ninsp, 0.); 811 Vector<float> rvecref(nchan_, 0.); 812 Vector<Complex> cvecref(nchan_/2+1, 0.); 813 uInt icol = 0; 814 unsigned int nreject = 0; 815 816 #ifdef KS_DEBUG 817 cout << "Starting initial FFT. The number of input spectra = " << ninsp << endl; 818 cout << "out matrix has ncolumn = " << outmat.ncolumn() << endl; 819 #endif 820 821 for (uInt isp = 0 ; isp < ninsp ; isp++) { 822 rvecref.reference( specmat.column(isp) ); 823 cvecref.reference( fftspmat.column(isp) ); 824 825 #ifdef KS_DEBUG 826 minMax(minval, maxval, rvecref); 827 cout << "Max/Min of inv FFTed Matrix = (max: " << maxval << ", min: " << minval << ")" << endl; 828 #endif 829 830 fftsf.fft0(cvecref, rvecref, true); 831 832 #ifdef KS_DEBUG 833 double maxr=cvecref[0].real(), minr=cvecref[0].real(), 834 maxi=cvecref[0].imag(), mini=cvecref[0].imag(); 835 for (uInt i = 1; i<cvecref.size();i++){ 836 maxr = max(maxr, cvecref[i].real()); 837 maxi = max(maxi, cvecref[i].imag()); 838 minr = min(minr, cvecref[i].real()); 839 mini = min(mini, cvecref[i].imag()); 840 } 841 cout << "Max/Min of inv FFTed Matrix (size=" << cvecref.size() << ") = (max: " << maxr << " + " << maxi << "j , min: " << minr << " + " << mini << "j)" << endl; 842 #endif 843 } 844 845 //Liberate from reference 846 rvecref.unique(); 847 848 Vector<Complex> cspec(nchan_/2+1, 0.); 849 const double PI = 6.0 * asin(0.5); 850 const double nchani = 1. / (float) nchan_; 851 const Complex trans(0., 1.); 852 #ifdef KS_DEBUG 853 cout << "starting actual deconvolution" << endl; 854 #endif 855 for (uInt j = 0 ; j < ninsp ; j++) { 856 for (uInt k = j+1 ; k < ninsp ; k++) { 857 const double dx = (shiftvec[k] - shiftvec[j]) * 2. * PI * nchani; 858 859 #ifdef KS_DEBUG 860 cout << "icol = " << icol << endl; 861 #endif 862 863 for (uInt ichan = 0 ; ichan < cspec.size() ; ichan++){ 864 cspec[ichan] = ( fftspmat(ichan, j) + fftspmat(ichan, k) )*0.5; 865 double phase = dx*ichan; 866 if ( fabs( sin(phase) ) > threshold){ 867 cspec[ichan] += ( fftspmat(ichan, j) - fftspmat(ichan, k) ) * 0.5 868 * trans * sin(phase) / ( 1. - cos(phase) ); 869 } else { 870 nreject++; 871 } 872 } // end of channel loop 873 874 #ifdef KS_DEBUG 875 cout << "done calculation of cspec" << endl; 876 #endif 877 878 Vector<Float> rspec; 879 rspec.reference( outmat.column(icol) ); 880 881 #ifdef KS_DEBUG 882 cout << "Starting inverse FFT. icol = " << icol << endl; 883 //cout << "- size of complex vector = " << cspec.size() << endl; 884 double maxr=cspec[0].real(), minr=cspec[0].real(), 885 maxi=cspec[0].imag(), mini=cspec[0].imag(); 886 for (uInt i = 1; i<cspec.size();i++){ 887 maxr = max(maxr, cspec[i].real()); 888 maxi = max(maxi, cspec[i].imag()); 889 minr = min(minr, cspec[i].real()); 890 mini = min(mini, cspec[i].imag()); 891 } 892 cout << "Max/Min of conv vector (size=" << cspec.size() << ") = (max: " << maxr << " + " << maxi << "j , min: " << minr << " + " << mini << "j)" << endl; 893 #endif 894 895 fftsi.fft0(rspec, cspec, false); 896 897 #ifdef KS_DEBUG 898 //cout << "- size of inversed real vector = " << rspec.size() << endl; 899 minMax(minval, maxval, rspec); 900 cout << "Max/Min of inv FFTed Vector (size=" << rspec.size() << ") = (max: " << maxval << ", min: " << minval << ")" << endl; 901 //cout << "Done inverse FFT. icol = " << icol << endl; 902 #endif 903 904 icol++; 905 } 906 } 907 908 #ifdef KS_DEBUG 909 minMax(minval, maxval, outmat); 910 cout << "Max/Min of inv FFTed Matrix = (max: " << maxval << ", min: " << minval << ")" << endl; 911 #endif 912 913 os << "Threshold = " << threshold << ", Rejected channels = " << nreject << endl; 914 }; 915 916 //////////////////////////////////////////////////////////////////// 917 // void STSideBandSep::cpprfft(std::vector<float> invec) 918 // { 919 // cout << "c++ method cpprfft" << endl; 920 // const unsigned int len = invec.size(); 921 // Vector<Complex> carr(len/2+1, 0.); 922 // Vector<float> inarr = Vector<float>(invec); 923 // Vector <float> outarr(len, 0.); 924 // FFTServer<Float, Complex> fftsf, fftsi; 925 // fftsf.resize(IPosition(1, len), FFTEnums::REALTOCOMPLEX); 926 // fftsi.resize(IPosition(1, invec.size()), FFTEnums::COMPLEXTOREAL); 927 // cout << "Input float array (length = " << len << "):" << endl; 928 // for (uInt i = 0 ; i < len ; i++){ 929 // cout << (i == 0 ? "( " : " ") << inarr[i] << (i == len-1 ? ")" : ","); 930 // } 931 // cout << endl; 932 // cout << "R->C transform" << endl; 933 // fftsf.fft0(carr, inarr, true); 934 // cout << "FFTed complex array (length = " << carr.size() << "):" << endl; 935 // for (uInt i = 0 ; i < carr.size() ; i++){ 936 // cout << (i == 0 ? "( " : " ") << carr[i] << ( (i == carr.size()-1) ? ")" : "," ); 937 // } 938 // cout << endl; 939 // cout << "C->R transform" << endl; 940 // fftsi.fft0(outarr, carr, false); 941 // cout << "invFFTed float array (length = " << outarr.size() << "):" << endl; 942 // for (uInt i = 0 ; i < outarr.size() ; i++){ 943 // cout << (i == 0 ? "( " : " ") << outarr[i] << ( (i == outarr.size()-1) ? ")" : "," ); 944 // } 945 // cout << endl; 946 // }; 947 //////////////////////////////////////////////////////////////////// 948 949 950 void STSideBandSep::aggregateMat(Matrix<float> &inmat, 951 vector<float> &outvec) 952 { 953 LogIO os(LogOrigin("STSideBandSep","aggregateMat()", WHERE)); 954 if (inmat.nrow() != nchan_) 955 throw(AipsError("Internal error. The row numbers of input matrix differs from nchan_")); 956 // if (outvec.size() != nchan_) 957 // throw(AipsError("Internal error. The size of output vector should be equal to nchan_")); 958 959 os << "Averaging " << inmat.ncolumn() << " spectra in the input matrix." 960 << LogIO::POST; 961 962 const uInt nspec = inmat.ncolumn(); 963 const double scale = 1./( (double) nspec ); 964 // initialize values with 0 965 outvec.assign(nchan_, 0); 966 for (uInt isp = 0 ; isp < nspec ; isp++) { 967 for (uInt ich = 0 ; ich < nchan_ ; ich++) { 968 outvec[ich] += inmat(ich, isp); 969 } 970 } 971 972 vector<float>::iterator iter; 973 for (iter = outvec.begin(); iter != outvec.end(); iter++){ 974 *iter *= scale; 975 } 976 }; 977 978 void STSideBandSep::subtractFromOther(const Matrix<float> &shiftmat, 979 const vector<float> &invec, 980 const vector<double> &shift, 981 vector<float> &outvec) 982 { 983 LogIO os(LogOrigin("STSideBandSep","subtractFromOther()", WHERE)); 984 if (shiftmat.nrow() != nchan_) 985 throw(AipsError("Internal error. The row numbers of input matrix differs from nchan_")); 986 if (invec.size() != nchan_) 987 throw(AipsError("Internal error. The length of input vector should be nchan_")); 988 if (shift.size() != shiftmat.ncolumn()) 989 throw(AipsError("Internal error. The column numbers of input matrix != the number of elements in shift")); 990 991 const uInt nspec = shiftmat.ncolumn(); 992 Vector<float> subsp(nchan_, 0.), shiftsub; 993 Matrix<float> submat(nchan_, nspec, 0.); 994 vector<float>::iterator iter; 995 for (uInt isp = 0 ; isp < nspec ; isp++) { 996 for (uInt ich = 0; ich < nchan_ ; ich++) { 997 subsp(ich) = shiftmat(ich, isp) - invec[ich]; 998 } 999 shiftsub.reference(submat.column(isp)); 1000 shiftSpectrum(subsp, shift[isp], shiftsub); 1001 } 1002 1003 aggregateMat(submat, outvec); 1004 }; 1005 1006 1007 void STSideBandSep::setLO1(const string lo1, const string frame, 249 1008 const double reftime, const string refdir) 250 1009 { 251 lo1Freq_ = lo1; 1010 Quantum<Double> qfreq; 1011 readQuantity(qfreq, String(lo1)); 1012 lo1Freq_ = qfreq.getValue("Hz"); 252 1013 MFrequency::getType(loFrame_, frame); 253 1014 loTime_ = reftime; … … 297 1058 }; 298 1059 299 ///// TEMPORAL FUNCTION!!! ///// 300 void STSideBandSep::setScanTb0(const ScantableWrapper &s){ 301 st0_ = s.getCP(); 302 }; 303 //////////////////////////////// 304 305 void STSideBandSep::solveImageFreqency() 306 { 307 #ifdef KS_DEBUG 308 cout << "STSideBandSep::solveImageFrequency" << endl; 309 #endif 1060 1061 void STSideBandSep::solveImageFrequency() 1062 { 310 1063 LogIO os(LogOrigin("STSideBandSep","solveImageFreqency()", WHERE)); 311 1064 os << "Start calculating frequencies of image side band" << LogIO::POST; … … 383 1136 sigrefval = toloframe(Quantum<Double>(refval, "Hz")).get("Hz").getValue(); 384 1137 385 386 1138 // Check for the availability of LO1 387 1139 if (lo1Freq_ > 0.) { … … 402 1154 // Try getting ASDM name from scantable header 403 1155 os << "Try getting information from scantable header" << LogIO::POST; 404 if (!getLo1FromScanTab( st0_, sigrefval, refpix, increment, nChan)) {1156 if (!getLo1FromScanTab(tableList_[0], sigrefval, refpix, increment, nChan)) { 405 1157 //throw AipsError("Failed to get LO1 frequency from asis table"); 406 1158 os << LogIO::WARN << "Failed to get LO1 frequency using information in scantable." << LogIO::POST; … … 410 1162 } 411 1163 } 1164 412 1165 // LO1 should now be ready. 413 1166 if (lo1Freq_ < 0.) 414 1167 throw(AipsError("Got negative LO1 Frequency")); 415 1168 416 // Print summary (LO1) 417 Vector<Double> dirvec = md.getAngle(Unit(String("rad"))).getValue(); 418 os << "[LO1 settings]" << LogIO::POST; 419 os << "- Frequency: " << lo1Freq_ << " [Hz] (" 420 << MFrequency::showType(loFrame_) << ")" << LogIO::POST; 421 os << "- Reference time: " << me.get(Unit(String("d"))).getValue() 422 << " [day]" << LogIO::POST; 423 os << "- Reference direction: [" << dirvec[0] << ", " << dirvec[1] 424 << "] (" << md.getRefString() << ") " << LogIO::POST; 425 426 //Print summary (signal) 427 os << "[Signal side band]" << LogIO::POST; 428 os << "- IFNO: " << imgTab_p->getIF(0) << " (FREQ_ID = " << freqid << ")" 429 << LogIO::POST; 430 os << "- Reference value: " << refval << " [Hz] (" 431 << MFrequency::showType(tabframe) << ") = " 432 << sigrefval << " [Hz] (" << MFrequency::showType(loFrame_) 433 << ")" << LogIO::POST; 434 os << "- Reference pixel: " << refpix << LogIO::POST; 435 os << "- Increment: " << increment << " [Hz]" << LogIO::POST; 1169 // Print summary 1170 { 1171 // LO1 1172 Vector<Double> dirvec = md.getAngle(Unit(String("rad"))).getValue(); 1173 os << "[LO1 settings]" << LogIO::POST; 1174 os << "- Frequency: " << lo1Freq_ << " [Hz] (" 1175 << MFrequency::showType(loFrame_) << ")" << LogIO::POST; 1176 os << "- Reference time: " << me.get(Unit(String("d"))).getValue() 1177 << " [day]" << LogIO::POST; 1178 os << "- Reference direction: [" << dirvec[0] << ", " << dirvec[1] 1179 << "] (" << md.getRefString() << ") " << LogIO::POST; 1180 1181 // signal sideband 1182 os << "[Signal side band]" << LogIO::POST; 1183 os << "- IFNO: " << imgTab_p->getIF(0) << " (FREQ_ID = " << freqid << ")" 1184 << LogIO::POST; 1185 os << "- Reference value: " << refval << " [Hz] (" 1186 << MFrequency::showType(tabframe) << ") = " 1187 << sigrefval << " [Hz] (" << MFrequency::showType(loFrame_) 1188 << ")" << LogIO::POST; 1189 os << "- Reference pixel: " << refpix << LogIO::POST; 1190 os << "- Increment: " << increment << " [Hz]" << LogIO::POST; 1191 } 436 1192 437 1193 // Calculate image band incr and refval in loFrame_ … … 447 1203 freqIdVec = fIDnew; 448 1204 449 // Print summary (Image side band) 450 os << "[Image side band]" << LogIO::POST; 451 os << "- IFNO: " << imgTab_p->getIF(0) << " (FREQ_ID = " << freqIdVec(0) 452 << ")" << LogIO::POST; 453 os << "- Reference value: " << imgrefval << " [Hz] (" 454 << MFrequency::showType(tabframe) << ") = " 455 << imgrefval_tab << " [Hz] (" << MFrequency::showType(loFrame_) 456 << ")" << LogIO::POST; 457 os << "- Reference pixel: " << refpix << LogIO::POST; 458 os << "- Increment: " << imgincr << " [Hz]" << LogIO::POST; 1205 // Print summary (Image sideband) 1206 { 1207 os << "[Image side band]" << LogIO::POST; 1208 os << "- IFNO: " << imgTab_p->getIF(0) << " (FREQ_ID = " << freqIdVec(0) 1209 << ")" << LogIO::POST; 1210 os << "- Reference value: " << imgrefval << " [Hz] (" 1211 << MFrequency::showType(tabframe) << ") = " 1212 << imgrefval_tab << " [Hz] (" << MFrequency::showType(loFrame_) 1213 << ")" << LogIO::POST; 1214 os << "- Reference pixel: " << refpix << LogIO::POST; 1215 os << "- Increment: " << imgincr << " [Hz]" << LogIO::POST; 1216 } 459 1217 }; 460 1218 … … 619 1377 }; 620 1378 621 // String STSideBandSep::()622 // {623 624 // };625 626 1379 } //namespace asap -
trunk/src/STSideBandSep.h
r2726 r2794 22 22 #include <coordinates/Coordinates/DirectionCoordinate.h> 23 23 #include <coordinates/Coordinates/SpectralCoordinate.h> 24 #include <scimath/Mathematics/FFTServer.h> 24 25 // asap 25 26 #include "ScantableWrapper.h" … … 40 41 explicit STSideBandSep(const vector<ScantableWrapper> &tables); 41 42 virtual ~STSideBandSep(); 43 44 ///////////// temp functions ////////////////////// 45 //void cpprfft(std::vector<float> invec); 46 ////////////////////////////////////////////////// 47 48 /** 49 * Separate side bands 50 **/ 51 void separate(string outname); 42 52 43 53 /** … … 83 93 * Set additional information to fill frequencies of image sideband 84 94 **/ 85 void setLO1(const doublelo1, const string frame="TOPO",95 void setLO1(const string lo1, const string frame="TOPO", 86 96 const double reftime=-1, string refdir=""); 87 97 void setLO1Root(const string name); 88 89 /**90 * Actual calculation of frequencies of image sideband91 **/92 void solveImageFreqency();93 98 94 99 private: … … 99 104 /** Return if the path exists (optionally, check file type) **/ 100 105 Bool checkFile(const string name, string type=""); 106 107 /** **/ 108 unsigned int setupShift(); 109 bool getFreqInfo(const CountedPtr<Scantable> &stab, const unsigned int &ifno, 110 double &freq0, double &incr, unsigned int &nchan); 111 112 /** Grid scantable **/ 113 ScantableWrapper gridTable(); 114 void mapExtent(vector< CountedPtr<Scantable> > &tablist, 115 Double &xmin, Double &xmax, 116 Double &ymin, Double &ymax); 117 118 /** 119 * Actual calculation of frequencies of image sideband 120 **/ 121 void solveImageFrequency(); 101 122 102 123 /** … … 112 133 const double refval, const double refpix, 113 134 const double increment, const int nChan); 135 bool getSpectraToSolve(const int polId, const int beamId, 136 const double dirX, const double dirY, 137 Matrix<float> &specmat, vector<uInt> &tabIdvec); 138 139 vector<float> solve(const Matrix<float> &specmat, 140 const vector<uInt> &tabIdvec, 141 const bool signal = true); 142 143 void shiftSpectrum(const Vector<float> &invec, double shift, 144 Vector<float> &outvec); 145 146 void deconvolve(Matrix<float> &specmat, const vector<double> shiftvec, 147 const double threshold, Matrix<float> &outmat); 148 149 void aggregateMat(Matrix<float> &inmat, vector<float> &outvec); 150 151 void subtractFromOther(const Matrix<float> &shiftmat, 152 const vector<float> &invec, 153 const vector<double> &shift, 154 vector<float> &outvec); 155 156 114 157 115 158 /** Member variables **/ … … 119 162 unsigned int ntable_; 120 163 // frequency and direction setup to select data. 121 unsignedint sigIfno_;164 int sigIfno_; 122 165 Quantum<Double> ftol_; 123 166 MFrequency::Types solFrame_; … … 130 173 double rejlimit_; 131 174 // LO1 132 double lo1Freq_; 175 double lo1Freq_; // in Hz 133 176 MFrequency::Types loFrame_; 134 177 double loTime_; … … 136 179 string asdmName_, asisName_; 137 180 181 //CountedPtr<Scantable> imgTab_p, sigTab_p; 138 182 CountedPtr<Scantable> imgTab_p, sigTab_p; 183 Table::TableType tp_; 184 FFTServer<Float, Complex> fftsf, fftsi; 139 185 // TEMPORAL member 140 186 CountedPtr<Scantable> st0_; -
trunk/src/python_STSideBandSep.cpp
r2726 r2794 34 34 boost::python::arg("refdir")="") ) 35 35 .def( "set_lo1root", &STSideBandSep::setLO1Root ) 36 // temporal methods 37 .def( "set_imgtable", &STSideBandSep::setImageTable ) 38 .def( "solve_imgfreq", &STSideBandSep::solveImageFreqency ) 39 .def( "_get_asistb_from_scantb", &STSideBandSep::setScanTb0 ) 36 .def( "separate", &STSideBandSep::separate ) 37 //// temporal methods 38 //.def( "_cpprfft", &STSideBandSep::cpprfft ) 39 //.def( "solve_imgfreq", &STSideBandSep::solveImageFreqency ) 40 //.def( "_get_asistb_from_scantb", &STSideBandSep::setScanTb0 ) 40 41 ; 41 42 };
Note:
See TracChangeset
for help on using the changeset viewer.