[2647] | 1 | import os, shutil
|
---|
| 2 | import numpy
|
---|
| 3 | import numpy.fft as FFT
|
---|
| 4 | import math
|
---|
| 5 |
|
---|
| 6 | from asap.scantable import scantable
|
---|
| 7 | from asap.parameters import rcParams
|
---|
| 8 | from asap.logging import asaplog, asaplog_post_dec
|
---|
| 9 | from asap.selector import selector
|
---|
| 10 | from asap.asapgrid import asapgrid2
|
---|
| 11 | #from asap._asap import sidebandsep
|
---|
| 12 |
|
---|
| 13 | class sbseparator:
|
---|
| 14 | """
|
---|
[2649] | 15 | The sbseparator class is defined to separate SIGNAL and IMAGE
|
---|
| 16 | sideband spectra observed by frequency-switching technique.
|
---|
| 17 | It also helps supressing emmission of IMAGE sideband.
|
---|
| 18 | *** WARNING *** THIS MODULE IS EXPERIMENTAL
|
---|
| 19 | Known issues:
|
---|
| 20 | - Frequencies of IMAGE sideband cannot be reconstructed from
|
---|
| 21 | information in scantable in sideband sparation. Frequency
|
---|
| 22 | setting of SIGNAL sideband is stored in output table for now.
|
---|
| 23 | - Flag information (stored in FLAGTRA) is ignored.
|
---|
[2647] | 24 |
|
---|
| 25 | Example:
|
---|
| 26 | # Create sideband separator instance whith 3 input data
|
---|
| 27 | sbsep = sbseparator(['test1.asap', 'test2.asap', 'test3.asap'])
|
---|
| 28 | # Set reference IFNO and tolerance to select data
|
---|
| 29 | sbsep.set_frequency(5, 30, frame='TOPO')
|
---|
| 30 | # Set direction tolerance to select data in unit of radian
|
---|
| 31 | sbsep.set_dirtol(1.e-5)
|
---|
| 32 | # Set rejection limit of solution
|
---|
| 33 | sbsep.set_limit(0.2)
|
---|
| 34 | # Solve image sideband as well
|
---|
| 35 | sbsep.set_both(True)
|
---|
| 36 | # Invoke sideband separation
|
---|
| 37 | sbsep.separate('testout.asap', overwrite = True)
|
---|
| 38 | """
|
---|
| 39 | def __init__(self, infiles):
|
---|
| 40 | self.intables = None
|
---|
| 41 | self.signalShift = []
|
---|
| 42 | self.imageShift = []
|
---|
| 43 | self.dsbmode = True
|
---|
| 44 | self.getboth = False
|
---|
| 45 | self.rejlimit = 0.2
|
---|
| 46 | self.baseif = -1
|
---|
| 47 | self.freqtol = 10.
|
---|
| 48 | self.freqframe = ""
|
---|
| 49 | self.solveother = False
|
---|
| 50 | self.dirtol = [1.e-5, 1.e-5] # direction tolerance in rad (2 arcsec)
|
---|
| 51 |
|
---|
| 52 | self.tables = []
|
---|
| 53 | self.nshift = -1
|
---|
| 54 | self.nchan = -1
|
---|
| 55 |
|
---|
| 56 | self.set_data(infiles)
|
---|
| 57 |
|
---|
| 58 | #self.separator = sidebandsep()
|
---|
| 59 |
|
---|
| 60 | @asaplog_post_dec
|
---|
| 61 | def set_data(self, infiles):
|
---|
| 62 | """
|
---|
| 63 | Set data to be processed.
|
---|
| 64 |
|
---|
| 65 | infiles : a list of filenames or scantables
|
---|
| 66 | """
|
---|
| 67 | if not (type(infiles) in (list, tuple, numpy.ndarray)):
|
---|
| 68 | infiles = [infiles]
|
---|
| 69 | if isinstance(infiles[0], scantable):
|
---|
| 70 | # a list of scantable
|
---|
| 71 | for stab in infiles:
|
---|
| 72 | if not isinstance(stab, scantable):
|
---|
| 73 | asaplog.post()
|
---|
| 74 | raise TypeError, "Input data is not a list of scantables."
|
---|
| 75 |
|
---|
| 76 | #self.separator._setdata(infiles)
|
---|
| 77 | self._reset_data()
|
---|
| 78 | self.intables = infiles
|
---|
| 79 | else:
|
---|
| 80 | # a list of filenames
|
---|
| 81 | for name in infiles:
|
---|
| 82 | if not os.path.exists(name):
|
---|
| 83 | asaplog.post()
|
---|
| 84 | raise ValueError, "Could not find input file '%s'" % name
|
---|
| 85 |
|
---|
| 86 | #self.separator._setdataname(infiles)
|
---|
| 87 | self._reset_data()
|
---|
| 88 | self.intables = infiles
|
---|
| 89 |
|
---|
| 90 | asaplog.push("%d files are set to process" % len(self.intables))
|
---|
| 91 |
|
---|
| 92 |
|
---|
| 93 | def _reset_data(self):
|
---|
| 94 | del self.intables
|
---|
| 95 | self.intables = None
|
---|
| 96 | self.signalShift = []
|
---|
| 97 | #self.imageShift = []
|
---|
| 98 | self.tables = []
|
---|
| 99 | self.nshift = -1
|
---|
| 100 | self.nchan = -1
|
---|
| 101 |
|
---|
| 102 | @asaplog_post_dec
|
---|
| 103 | def set_frequency(self, baseif, freqtol, frame=""):
|
---|
| 104 | """
|
---|
| 105 | Set IFNO and frequency tolerance to select data to process.
|
---|
| 106 |
|
---|
| 107 | Parameters:
|
---|
| 108 | - reference IFNO to process in the first table in the list
|
---|
| 109 | - frequency tolerance from reference IF to select data
|
---|
| 110 | frame : frequency frame to select IF
|
---|
| 111 | """
|
---|
| 112 | self._reset_if()
|
---|
| 113 | self.baseif = baseif
|
---|
| 114 | if isinstance(freqtol,dict) and freqtol["unit"] == "Hz":
|
---|
| 115 | if freqtol['value'] > 0.:
|
---|
| 116 | self.freqtol = freqtol
|
---|
| 117 | else:
|
---|
| 118 | asaplog.post()
|
---|
| 119 | asaplog.push("Frequency tolerance should be positive value.")
|
---|
| 120 | asaplog.post("ERROR")
|
---|
| 121 | return
|
---|
| 122 | else:
|
---|
| 123 | # torelance in channel unit
|
---|
| 124 | if freqtol > 0:
|
---|
| 125 | self.freqtol = float(freqtol)
|
---|
| 126 | else:
|
---|
| 127 | asaplog.post()
|
---|
| 128 | asaplog.push("Frequency tolerance should be positive value.")
|
---|
| 129 | asaplog.post("ERROR")
|
---|
| 130 | return
|
---|
| 131 | self.freqframe = frame
|
---|
| 132 |
|
---|
| 133 | def _reset_if(self):
|
---|
| 134 | self.baseif = -1
|
---|
| 135 | self.freqtol = 0
|
---|
| 136 | self.freqframe = ""
|
---|
| 137 | self.signalShift = []
|
---|
| 138 | #self.imageShift = []
|
---|
| 139 | self.tables = []
|
---|
| 140 | self.nshift = 0
|
---|
| 141 | self.nchan = -1
|
---|
| 142 |
|
---|
| 143 | @asaplog_post_dec
|
---|
| 144 | def set_dirtol(self, dirtol=[1.e-5,1.e-5]):
|
---|
| 145 | """
|
---|
| 146 | Set tolerance of direction to select data
|
---|
| 147 | """
|
---|
| 148 | # direction tolerance in rad
|
---|
| 149 | if not (type(dirtol) in [list, tuple, numpy.ndarray]):
|
---|
| 150 | dirtol = [dirtol, dirtol]
|
---|
| 151 | if len(dirtol) == 1:
|
---|
| 152 | dirtol = [dirtol[0], dirtol[0]]
|
---|
| 153 | if len(dirtol) > 1:
|
---|
| 154 | self.dirtol = dirtol[0:2]
|
---|
| 155 | else:
|
---|
| 156 | asaplog.post()
|
---|
| 157 | asaplog.push("Invalid direction tolerance. Should be a list of float in unit radian")
|
---|
| 158 | asaplog.post("ERROR")
|
---|
| 159 | return
|
---|
| 160 | asaplog.post("Set direction tolerance [%f, %f] (rad)" % \
|
---|
| 161 | (self.dirtol[0], self.dirtol[1]))
|
---|
| 162 |
|
---|
| 163 | @asaplog_post_dec
|
---|
| 164 | def set_shift(self, mode="DSB", imageshift=None):
|
---|
| 165 | """
|
---|
| 166 | Set shift mode and channel shift of image band.
|
---|
| 167 |
|
---|
| 168 | mode : shift mode ['DSB'|'SSB']
|
---|
| 169 | When mode='DSB', imageshift is assumed to be equal
|
---|
| 170 | to the shift of signal sideband but in opposite direction.
|
---|
| 171 | imageshift : a list of number of channel shift in image band of
|
---|
| 172 | each scantable. valid only mode='SSB'
|
---|
| 173 | """
|
---|
| 174 | if mode.upper().startswith("S"):
|
---|
| 175 | if not imageshift:
|
---|
| 176 | raise ValueError, "Need to set shift value of image sideband"
|
---|
| 177 | self.dsbmode = False
|
---|
| 178 | self.imageShift = imageshift
|
---|
| 179 | asaplog.push("Image sideband shift is set manually: %s" % str(self.imageShift))
|
---|
| 180 | else:
|
---|
| 181 | # DSB mode
|
---|
| 182 | self.dsbmode = True
|
---|
| 183 | self.imageShift = []
|
---|
| 184 |
|
---|
| 185 | @asaplog_post_dec
|
---|
| 186 | def set_both(self, flag=False):
|
---|
| 187 | """
|
---|
| 188 | Resolve both image and signal sideband when True.
|
---|
| 189 | """
|
---|
| 190 | self.getboth = flag
|
---|
| 191 | if self.getboth:
|
---|
| 192 | asaplog.push("Both signal and image sidebands are solved and output as separate tables.")
|
---|
| 193 | else:
|
---|
| 194 | asaplog.push("Only signal sideband is solved and output as an table.")
|
---|
| 195 |
|
---|
| 196 | @asaplog_post_dec
|
---|
| 197 | def set_limit(self, threshold=0.2):
|
---|
| 198 | """
|
---|
| 199 | Set rejection limit of solution.
|
---|
| 200 | """
|
---|
| 201 | #self.separator._setlimit(abs(threshold))
|
---|
| 202 | self.rejlimit = threshold
|
---|
| 203 | asaplog.push("The threshold of rejection is set to %f" % self.rejlimit)
|
---|
| 204 |
|
---|
| 205 |
|
---|
| 206 | @asaplog_post_dec
|
---|
| 207 | def set_solve_other(self, flag=False):
|
---|
| 208 | """
|
---|
| 209 | Calculate spectra by subtracting the solution of the other sideband
|
---|
| 210 | when True.
|
---|
| 211 | """
|
---|
| 212 | self.solveother = flag
|
---|
| 213 | if flag:
|
---|
| 214 | asaplog.push("Expert mode: solution are obtained by subtraction of the other sideband.")
|
---|
| 215 |
|
---|
| 216 |
|
---|
| 217 | @asaplog_post_dec
|
---|
| 218 | def separate(self, outname="", overwrite=False):
|
---|
| 219 | """
|
---|
| 220 | Invoke sideband separation.
|
---|
| 221 |
|
---|
| 222 | outname : a name of output scantable
|
---|
| 223 | overwrite : overwrite existing table
|
---|
| 224 | """
|
---|
| 225 | # List up valid scantables and IFNOs to convolve.
|
---|
| 226 | #self.separator._separate()
|
---|
| 227 | self._setup_shift()
|
---|
| 228 | #self._preprocess_tables()
|
---|
| 229 |
|
---|
| 230 | nshift = len(self.tables)
|
---|
| 231 | signaltab = self._grid_outtable(self.tables[0].copy())
|
---|
| 232 | if self.getboth:
|
---|
| 233 | imagetab = signaltab.copy()
|
---|
| 234 |
|
---|
| 235 | rejrow = []
|
---|
| 236 | for irow in xrange(signaltab.nrow()):
|
---|
| 237 | currpol = signaltab.getpol(irow)
|
---|
| 238 | currbeam = signaltab.getbeam(irow)
|
---|
| 239 | currdir = signaltab.get_directionval(irow)
|
---|
| 240 | spec_array, tabidx = self._get_specarray(polid=currpol,\
|
---|
| 241 | beamid=currbeam,\
|
---|
| 242 | dir=currdir)
|
---|
| 243 | #if not spec_array:
|
---|
| 244 | if len(tabidx) == 0:
|
---|
| 245 | asaplog.post()
|
---|
| 246 | asaplog.push("skipping row %d" % irow)
|
---|
| 247 | rejrow.append(irow)
|
---|
| 248 | continue
|
---|
| 249 | signal = self._solve_signal(spec_array, tabidx)
|
---|
| 250 | signaltab.set_spectrum(signal, irow)
|
---|
| 251 | if self.getboth:
|
---|
| 252 | image = self._solve_image(spec_array, tabidx)
|
---|
| 253 | imagetab.set_spectrum(image, irow)
|
---|
| 254 |
|
---|
| 255 | # TODO: Need to remove rejrow form scantables here
|
---|
| 256 | signaltab.flag_row(rejrow)
|
---|
| 257 | if self.getboth:
|
---|
| 258 | imagetab.flag_row(rejrow)
|
---|
| 259 |
|
---|
| 260 | if outname == "":
|
---|
| 261 | outname = "sbsepareted.asap"
|
---|
| 262 | signalname = outname + ".signalband"
|
---|
| 263 | if os.path.exists(signalname):
|
---|
| 264 | if not overwrite:
|
---|
| 265 | raise Exception, "Output file '%s' exists." % signalname
|
---|
| 266 | else:
|
---|
| 267 | shutil.rmtree(signalname)
|
---|
| 268 | signaltab.save(signalname)
|
---|
| 269 | if self.getboth:
|
---|
[2649] | 270 | # Warnings
|
---|
| 271 | asaplog.post()
|
---|
| 272 | asaplog.push("Saving IMAGE sideband.")
|
---|
| 273 | asaplog.push("Note, frequency information of IMAGE sideband cannot be properly filled so far. (future development)")
|
---|
| 274 | asaplog.push("Storing frequency setting of SIGNAL sideband in FREQUENCIES table for now.")
|
---|
| 275 | asaplog.post("WARN")
|
---|
| 276 |
|
---|
[2647] | 277 | imagename = outname + ".imageband"
|
---|
| 278 | if os.path.exists(imagename):
|
---|
| 279 | if not overwrite:
|
---|
| 280 | raise Exception, "Output file '%s' exists." % imagename
|
---|
| 281 | else:
|
---|
| 282 | shutil.rmtree(imagename)
|
---|
| 283 | imagetab.save(imagename)
|
---|
| 284 |
|
---|
| 285 |
|
---|
| 286 | def _solve_signal(self, data, tabidx=None):
|
---|
| 287 | if not tabidx:
|
---|
| 288 | tabidx = range(len(data))
|
---|
| 289 |
|
---|
| 290 | tempshift = []
|
---|
| 291 | dshift = []
|
---|
| 292 | if self.solveother:
|
---|
| 293 | for idx in tabidx:
|
---|
| 294 | tempshift.append(-self.imageShift[idx])
|
---|
| 295 | dshift.append(self.signalShift[idx] - self.imageShift[idx])
|
---|
| 296 | else:
|
---|
| 297 | for idx in tabidx:
|
---|
| 298 | tempshift.append(-self.signalShift[idx])
|
---|
| 299 | dshift.append(self.imageShift[idx] - self.signalShift[idx])
|
---|
| 300 |
|
---|
| 301 | shiftdata = numpy.zeros(data.shape, numpy.float)
|
---|
| 302 | for i in range(len(data)):
|
---|
| 303 | shiftdata[i] = self._shiftSpectrum(data[i], tempshift[i])
|
---|
| 304 | ifftdata = self._Deconvolution(shiftdata, dshift, self.rejlimit)
|
---|
| 305 | result_image = self._combineResult(ifftdata)
|
---|
| 306 | if not self.solveother:
|
---|
| 307 | return result_image
|
---|
| 308 | result_signal = self._subtractOtherSide(shiftdata, dshift, result_image)
|
---|
| 309 | return result_signal
|
---|
| 310 |
|
---|
| 311 |
|
---|
| 312 | def _solve_image(self, data, tabidx=None):
|
---|
| 313 | if not tabidx:
|
---|
| 314 | tabidx = range(len(data))
|
---|
| 315 |
|
---|
| 316 | tempshift = []
|
---|
| 317 | dshift = []
|
---|
| 318 | if self.solveother:
|
---|
| 319 | for idx in tabidx:
|
---|
| 320 | tempshift.append(-self.signalShift[idx])
|
---|
| 321 | dshift.append(self.imageShift[idx] - self.signalShift[idx])
|
---|
| 322 | else:
|
---|
| 323 | for idx in tabidx:
|
---|
| 324 | tempshift.append(-self.imageShift[idx])
|
---|
| 325 | dshift.append(self.signalShift[idx] - self.imageShift[idx])
|
---|
| 326 |
|
---|
| 327 | shiftdata = numpy.zeros(data.shape, numpy.float)
|
---|
| 328 | for i in range(len(data)):
|
---|
| 329 | shiftdata[i] = self._shiftSpectrum(data[i], tempshift[i])
|
---|
| 330 | ifftdata = self._Deconvolution(shiftdata, dshift, self.rejlimit)
|
---|
| 331 | result_image = self._combineResult(ifftdata)
|
---|
| 332 | if not self.solveother:
|
---|
| 333 | return result_image
|
---|
| 334 | result_signal = self._subtractOtherSide(shiftdata, dshift, result_image)
|
---|
| 335 | return result_signal
|
---|
| 336 |
|
---|
| 337 | @asaplog_post_dec
|
---|
| 338 | def _grid_outtable(self, table):
|
---|
| 339 | # Generate gridded table for output (Just to get rows)
|
---|
| 340 | gridder = asapgrid2(table)
|
---|
| 341 | gridder.setIF(self.baseif)
|
---|
| 342 |
|
---|
| 343 | cellx = str(self.dirtol[0])+"rad"
|
---|
| 344 | celly = str(self.dirtol[1])+"rad"
|
---|
| 345 | dirarr = numpy.array(table.get_directionval()).transpose()
|
---|
| 346 | mapx = dirarr[0].max() - dirarr[0].min()
|
---|
| 347 | mapy = dirarr[1].max() - dirarr[1].min()
|
---|
| 348 | nx = max(1, numpy.ceil(mapx/self.dirtol[0]))
|
---|
| 349 | ny = max(1, numpy.ceil(mapy/self.dirtol[0]))
|
---|
[2685] | 350 |
|
---|
[2647] | 351 | asaplog.push("Regrid output scantable with cell = [%s, %s]" % \
|
---|
| 352 | (cellx, celly))
|
---|
| 353 | gridder.defineImage(nx=nx, ny=ny, cellx=cellx, celly=celly)
|
---|
[2685] | 354 | gridder.setFunc(func='box', convsupport=1)
|
---|
[2647] | 355 | gridder.setWeight(weightType='uniform')
|
---|
| 356 | gridder.grid()
|
---|
| 357 | return gridder.getResult()
|
---|
| 358 |
|
---|
| 359 | @asaplog_post_dec
|
---|
| 360 | def _get_specarray(self, polid=None, beamid=None, dir=None):
|
---|
| 361 | ntable = len(self.tables)
|
---|
| 362 | spec_array = numpy.zeros((ntable, self.nchan), numpy.float)
|
---|
| 363 | nspec = 0
|
---|
| 364 | asaplog.push("Start data selection by POL=%d, BEAM=%d, direction=[%f, %f]" % (polid, beamid, dir[0], dir[1]))
|
---|
| 365 | tabidx = []
|
---|
| 366 | for itab in range(ntable):
|
---|
| 367 | tab = self.tables[itab]
|
---|
| 368 | # Select rows by POLNO and BEAMNO
|
---|
| 369 | try:
|
---|
| 370 | tab.set_selection(pols=[polid], beams=[beamid])
|
---|
| 371 | if tab.nrow() > 0: tabidx.append(itab)
|
---|
| 372 | except: # no selection
|
---|
| 373 | asaplog.post()
|
---|
| 374 | asaplog.push("table %d - No spectrum ....skipping the table" % (itab))
|
---|
| 375 | asaplog.post("WARN")
|
---|
| 376 | continue
|
---|
| 377 |
|
---|
| 378 | # Select rows by direction
|
---|
| 379 | spec = numpy.zeros(self.nchan, numpy.float)
|
---|
| 380 | selrow = []
|
---|
| 381 | for irow in range(tab.nrow()):
|
---|
| 382 | currdir = tab.get_directionval(irow)
|
---|
| 383 | if (abs(currdir[0] - dir[0]) > self.dirtol[0]) or \
|
---|
| 384 | (abs(currdir[1] - dir[1]) > self.dirtol[1]):
|
---|
| 385 | continue
|
---|
| 386 | selrow.append(irow)
|
---|
| 387 | if len(selrow) == 0:
|
---|
| 388 | asaplog.post()
|
---|
| 389 | asaplog.push("table %d - No spectrum ....skipping the table" % (itab))
|
---|
| 390 | asaplog.post("WARN")
|
---|
| 391 | continue
|
---|
| 392 | else:
|
---|
| 393 | seltab = tab.copy()
|
---|
| 394 | seltab.set_selection(selector(rows=selrow))
|
---|
| 395 |
|
---|
| 396 | if tab.nrow() > 1:
|
---|
| 397 | asaplog.push("table %d - More than a spectrum selected. averaging rows..." % (itab))
|
---|
| 398 | tab = seltab.average_time(scanav=False, weight="tintsys")
|
---|
| 399 | else:
|
---|
| 400 | tab = seltab
|
---|
| 401 |
|
---|
| 402 | spec_array[nspec] = tab._getspectrum()
|
---|
| 403 | nspec += 1
|
---|
| 404 |
|
---|
| 405 | if nspec != ntable:
|
---|
| 406 | asaplog.post()
|
---|
| 407 | #asaplog.push("Some tables has no spectrum with POL=%d BEAM=%d. averaging rows..." % (polid, beamid))
|
---|
| 408 | asaplog.push("Could not find corresponding rows in some tables.")
|
---|
| 409 | asaplog.push("Number of spectra selected = %d (table: %d)" % (nspec, ntable))
|
---|
| 410 | if nspec < 2:
|
---|
| 411 | asaplog.push("At least 2 spectra are necessary for convolution")
|
---|
| 412 | asaplog.post("ERROR")
|
---|
| 413 | return False, tabidx
|
---|
| 414 |
|
---|
| 415 | return spec_array[0:nspec], tabidx
|
---|
| 416 |
|
---|
| 417 |
|
---|
| 418 | @asaplog_post_dec
|
---|
| 419 | def _setup_shift(self):
|
---|
| 420 | ### define self.tables, self.signalShift, and self.imageShift
|
---|
| 421 | if not self.intables:
|
---|
| 422 | asaplog.post()
|
---|
| 423 | raise RuntimeError, "Input data is not defined."
|
---|
| 424 | #if self.baseif < 0:
|
---|
| 425 | # asaplog.post()
|
---|
| 426 | # raise RuntimeError, "Reference IFNO is not defined."
|
---|
| 427 |
|
---|
| 428 | byname = False
|
---|
| 429 | #if not self.intables:
|
---|
| 430 | if isinstance(self.intables[0], str):
|
---|
| 431 | # A list of file name is given
|
---|
| 432 | if not os.path.exists(self.intables[0]):
|
---|
| 433 | asaplog.post()
|
---|
| 434 | raise RuntimeError, "Could not find '%s'" % self.intables[0]
|
---|
| 435 |
|
---|
| 436 | stab = scantable(self.intables[0],average=False)
|
---|
| 437 | ntab = len(self.intables)
|
---|
| 438 | byname = True
|
---|
| 439 | else:
|
---|
| 440 | stab = self.intables[0]
|
---|
| 441 | ntab = len(self.intables)
|
---|
| 442 |
|
---|
| 443 | if len(stab.getbeamnos()) > 1:
|
---|
| 444 | asaplog.post()
|
---|
| 445 | asaplog.push("Mult-beam data is not supported by this module.")
|
---|
| 446 | asaplog.post("ERROR")
|
---|
| 447 | return
|
---|
| 448 |
|
---|
| 449 | valid_ifs = stab.getifnos()
|
---|
| 450 | if self.baseif < 0:
|
---|
| 451 | self.baseif = valid_ifs[0]
|
---|
| 452 | asaplog.post()
|
---|
| 453 | asaplog.push("IFNO is not selected. Using the first IF in the first scantable. Reference IFNO = %d" % (self.baseif))
|
---|
| 454 |
|
---|
| 455 | if not (self.baseif in valid_ifs):
|
---|
| 456 | asaplog.post()
|
---|
| 457 | errmsg = "IF%d does not exist in the first scantable" % \
|
---|
| 458 | self.baseif
|
---|
| 459 | raise RuntimeError, errmsg
|
---|
| 460 |
|
---|
| 461 | asaplog.push("Start selecting tables and IFNOs to solve.")
|
---|
| 462 | asaplog.push("Cheching frequency of the reference IF")
|
---|
| 463 | unit_org = stab.get_unit()
|
---|
| 464 | coord = stab._getcoordinfo()
|
---|
| 465 | frame_org = coord[1]
|
---|
| 466 | stab.set_unit("Hz")
|
---|
| 467 | if len(self.freqframe) > 0:
|
---|
| 468 | stab.set_freqframe(self.freqframe)
|
---|
| 469 | stab.set_selection(ifs=[self.baseif])
|
---|
| 470 | spx = stab._getabcissa()
|
---|
| 471 | stab.set_selection()
|
---|
| 472 | basech0 = spx[0]
|
---|
| 473 | baseinc = spx[1]-spx[0]
|
---|
| 474 | self.nchan = len(spx)
|
---|
| 475 | if isinstance(self.freqtol, float):
|
---|
| 476 | vftol = abs(baseinc * self.freqtol)
|
---|
| 477 | self.freqtol = dict(value=vftol, unit="Hz")
|
---|
| 478 | else:
|
---|
| 479 | vftol = abs(self.freqtol['value'])
|
---|
| 480 | inctol = abs(baseinc/float(self.nchan))
|
---|
| 481 | asaplog.push("Reference frequency setup (Table = 0, IFNO = %d): nchan = %d, chan0 = %f Hz, incr = %f Hz" % (self.baseif, self.nchan, basech0, baseinc))
|
---|
| 482 | asaplog.push("Allowed frequency tolerance = %f Hz ( %f channels)" % (vftol, vftol/baseinc))
|
---|
| 483 | poltype0 = stab.poltype()
|
---|
| 484 |
|
---|
| 485 | self.tables = []
|
---|
| 486 | self.signalShift = []
|
---|
| 487 | if self.dsbmode:
|
---|
| 488 | self.imageShift = []
|
---|
| 489 |
|
---|
| 490 | for itab in range(ntab):
|
---|
| 491 | asaplog.push("Table %d:" % itab)
|
---|
| 492 | tab_selected = False
|
---|
| 493 | if itab > 0:
|
---|
| 494 | if byname:
|
---|
| 495 | stab = scantable(self.intables[itab],average=False)
|
---|
| 496 | else:
|
---|
| 497 | stab = self.intables[itab]
|
---|
| 498 | unit_org = stab.get_unit()
|
---|
| 499 | coord = stab._getcoordinfo()
|
---|
| 500 | frame_org = coord[1]
|
---|
| 501 | stab.set_unit("Hz")
|
---|
| 502 | if len(self.freqframe) > 0:
|
---|
| 503 | stab.set_freqframe(self.freqframe)
|
---|
| 504 |
|
---|
| 505 | # Check POLTYPE should be equal to the first table.
|
---|
| 506 | if stab.poltype() != poltype0:
|
---|
| 507 | asaplog.post()
|
---|
| 508 | raise Exception, "POLTYPE should be equal to the first table."
|
---|
| 509 | # Multiple beam data may not handled properly
|
---|
| 510 | if len(stab.getbeamnos()) > 1:
|
---|
| 511 | asaplog.post()
|
---|
| 512 | asaplog.push("table contains multiple beams. It may not be handled properly.")
|
---|
| 513 | asaplog.push("WARN")
|
---|
[2685] | 514 |
|
---|
[2647] | 515 | for ifno in stab.getifnos():
|
---|
| 516 | stab.set_selection(ifs=[ifno])
|
---|
| 517 | spx = stab._getabcissa()
|
---|
| 518 | if (len(spx) != self.nchan) or \
|
---|
| 519 | (abs(spx[0]-basech0) > vftol) or \
|
---|
| 520 | (abs(spx[1]-spx[0]-baseinc) > inctol):
|
---|
| 521 | continue
|
---|
| 522 | tab_selected = True
|
---|
| 523 | seltab = stab.copy()
|
---|
| 524 | seltab.set_unit(unit_org)
|
---|
| 525 | seltab.set_freqframe(frame_org)
|
---|
| 526 | self.tables.append(seltab)
|
---|
| 527 | self.signalShift.append((spx[0]-basech0)/baseinc)
|
---|
| 528 | if self.dsbmode:
|
---|
| 529 | self.imageShift.append(-self.signalShift[-1])
|
---|
| 530 | asaplog.push("- IF%d selected: sideband shift = %16.12e channels" % (ifno, self.signalShift[-1]))
|
---|
| 531 | stab.set_selection()
|
---|
| 532 | stab.set_unit(unit_org)
|
---|
| 533 | stab.set_freqframe(frame_org)
|
---|
| 534 | if not tab_selected:
|
---|
| 535 | asaplog.post()
|
---|
| 536 | asaplog.push("No data selected in Table %d" % itab)
|
---|
| 537 | asaplog.post("WARN")
|
---|
| 538 |
|
---|
| 539 | asaplog.push("Total number of IFs selected = %d" % len(self.tables))
|
---|
| 540 | if len(self.tables) < 2:
|
---|
| 541 | asaplog.post()
|
---|
| 542 | raise RuntimeError, "At least 2 IFs are necessary for convolution!"
|
---|
| 543 |
|
---|
| 544 | if not self.dsbmode and len(self.imageShift) != len(self.signalShift):
|
---|
| 545 | asaplog.post()
|
---|
| 546 | errmsg = "User defined channel shift of image sideband has %d elements, while selected IFNOs are %d" % (len(self.imageShift), len(self.signalShift))
|
---|
| 547 | raise RuntimeError, errmsg
|
---|
| 548 |
|
---|
| 549 | self.signalShift = numpy.array(self.signalShift)
|
---|
| 550 | self.imageShift = numpy.array(self.imageShift)
|
---|
| 551 | self.nshift = len(self.tables)
|
---|
| 552 |
|
---|
| 553 | @asaplog_post_dec
|
---|
| 554 | def _preprocess_tables(self):
|
---|
| 555 | ### temporary method to preprocess data
|
---|
| 556 | ### Do time averaging for now.
|
---|
| 557 | for itab in range(len(self.tables)):
|
---|
| 558 | self.tables[itab] = self.tables[itab].average_time(scanav=False, weight="tintsys")
|
---|
| 559 |
|
---|
| 560 |
|
---|
| 561 | # def save(self, outfile, outform="ASAP", overwrite=False):
|
---|
| 562 | # if not overwrite and os.path.exists(outfile):
|
---|
| 563 | # raise RuntimeError, "Output file '%s' already exists" % outfile
|
---|
| 564 | #
|
---|
| 565 | # #self.separator._save(outfile, outform)
|
---|
| 566 |
|
---|
| 567 | # def done(self):
|
---|
| 568 | # self.close()
|
---|
| 569 |
|
---|
| 570 | # def close(self):
|
---|
| 571 | # pass
|
---|
| 572 | # #del self.separator
|
---|
| 573 |
|
---|
| 574 |
|
---|
| 575 |
|
---|
| 576 | ########################################################################
|
---|
| 577 | def _Deconvolution(self, data_array, shift, threshold=0.00000001):
|
---|
| 578 | FObs = []
|
---|
| 579 | Reject = 0
|
---|
| 580 | nshift, nchan = data_array.shape
|
---|
| 581 | nspec = nshift*(nshift-1)/2
|
---|
| 582 | ifftObs = numpy.zeros((nspec, nchan), numpy.float)
|
---|
| 583 | for i in range(nshift):
|
---|
| 584 | F = FFT.fft(data_array[i])
|
---|
| 585 | FObs.append(F)
|
---|
| 586 | z = 0
|
---|
| 587 | for i in range(nshift):
|
---|
| 588 | for j in range(i+1, nshift):
|
---|
| 589 | Fobs = (FObs[i]+FObs[j])/2.0
|
---|
| 590 | dX = (shift[j]-shift[i])*2.0*math.pi/float(self.nchan)
|
---|
| 591 | #print 'dX,i,j=',dX,i,j
|
---|
| 592 | for k in range(1,self.nchan):
|
---|
| 593 | if math.fabs(math.sin(dX*k)) > threshold:
|
---|
| 594 | Fobs[k] += ((FObs[i][k]-FObs[j][k])/2.0/(1.0-math.cos(dX*k))*math.sin(dX*k))*1.0j
|
---|
| 595 | else: Reject += 1
|
---|
| 596 | ifftObs[z] = FFT.ifft(Fobs)
|
---|
| 597 | z += 1
|
---|
| 598 | print 'Threshold=%s Reject=%d' % (threshold, Reject)
|
---|
| 599 | return ifftObs
|
---|
| 600 |
|
---|
| 601 | def _combineResult(self, ifftObs):
|
---|
| 602 | nspec = len(ifftObs)
|
---|
| 603 | sum = ifftObs[0]
|
---|
| 604 | for i in range(1,nspec):
|
---|
| 605 | sum += ifftObs[i]
|
---|
| 606 | return(sum/float(nspec))
|
---|
| 607 |
|
---|
| 608 | def _subtractOtherSide(self, data_array, shift, Data):
|
---|
| 609 | sum = numpy.zeros(len(Data), numpy.float)
|
---|
| 610 | numSP = len(data_array)
|
---|
| 611 | for i in range(numSP):
|
---|
| 612 | SPsub = data_array[i] - Data
|
---|
| 613 | sum += self._shiftSpectrum(SPsub, -shift[i])
|
---|
| 614 | return(sum/float(numSP))
|
---|
| 615 |
|
---|
| 616 | def _shiftSpectrum(self, data, Shift):
|
---|
| 617 | Out = numpy.zeros(self.nchan, numpy.float)
|
---|
| 618 | w2 = Shift % 1
|
---|
| 619 | w1 = 1.0 - w2
|
---|
| 620 | for i in range(self.nchan):
|
---|
| 621 | c1 = int((Shift + i) % self.nchan)
|
---|
| 622 | c2 = (c1 + 1) % self.nchan
|
---|
| 623 | Out[c1] += data[i] * w1
|
---|
| 624 | Out[c2] += data[i] * w2
|
---|
| 625 | return Out.copy()
|
---|