| 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 SBSeparator | 
|---|
| 12 |  | 
|---|
| 13 | class sbseparator: | 
|---|
| 14 | """ | 
|---|
| 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. | 
|---|
| 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 | #self.lo1 = 0. | 
|---|
| 52 |  | 
|---|
| 53 | self.tables = [] | 
|---|
| 54 | self.nshift = -1 | 
|---|
| 55 | self.nchan = -1 | 
|---|
| 56 |  | 
|---|
| 57 | self.set_data(infiles) | 
|---|
| 58 |  | 
|---|
| 59 | self._separator = SBSeparator() | 
|---|
| 60 |  | 
|---|
| 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 | 
|---|
| 104 | def set_frequency(self, baseif, freqtol, frame=""): | 
|---|
| 105 | """ | 
|---|
| 106 | Set IFNO and frequency tolerance to select data to process. | 
|---|
| 107 |  | 
|---|
| 108 | Parameters: | 
|---|
| 109 | - reference IFNO to process in the first table in the list | 
|---|
| 110 | - frequency tolerance from reference IF to select data | 
|---|
| 111 | frame  : frequency frame to select IF | 
|---|
| 112 | """ | 
|---|
| 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 | 
|---|
| 118 | else: | 
|---|
| 119 | asaplog.post() | 
|---|
| 120 | asaplog.push("Frequency tolerance should be positive value.") | 
|---|
| 121 | 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) | 
|---|
| 294 |  | 
|---|
| 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." % signalname | 
|---|
| 301 | else: | 
|---|
| 302 | shutil.rmtree(signalname) | 
|---|
| 303 | signaltab.save(signalname) | 
|---|
| 304 |  | 
|---|
| 305 | if self.getboth: | 
|---|
| 306 | # Warnings | 
|---|
| 307 | 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." % imagename | 
|---|
| 317 | else: | 
|---|
| 318 | shutil.rmtree(imagename) | 
|---|
| 319 | # Update frequency information | 
|---|
| 320 | 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_image | 
|---|
| 347 | result_signal = self._subtractOtherSide(shiftdata, dshift, result_image) | 
|---|
| 348 | return result_signal | 
|---|
| 349 |  | 
|---|
| 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_image | 
|---|
| 373 | result_signal = self._subtractOtherSide(shiftdata, dshift, result_image) | 
|---|
| 374 | return result_signal | 
|---|
| 375 |  | 
|---|
| 376 | @asaplog_post_dec | 
|---|
| 377 | 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_dec | 
|---|
| 400 | 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 = 0 | 
|---|
| 404 | 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 BEAMNO | 
|---|
| 409 | try: | 
|---|
| 410 | tab.set_selection(pols=[polid], beams=[beamid]) | 
|---|
| 411 | if tab.nrow() > 0: tabidx.append(itab) | 
|---|
| 412 | except: # no selection | 
|---|
| 413 | asaplog.post() | 
|---|
| 414 | asaplog.push("table %d - No spectrum ....skipping the table" % (itab)) | 
|---|
| 415 | asaplog.post("WARN") | 
|---|
| 416 | continue | 
|---|
| 417 |  | 
|---|
| 418 | # Select rows by direction | 
|---|
| 419 | 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 | continue | 
|---|
| 426 | 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 | continue | 
|---|
| 432 | else: | 
|---|
| 433 | seltab = tab.copy() | 
|---|
| 434 | seltab.set_selection(selector(rows=selrow)) | 
|---|
| 435 |  | 
|---|
| 436 | if tab.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 = seltab | 
|---|
| 441 |  | 
|---|
| 442 | spec_array[nspec] = tab._getspectrum() | 
|---|
| 443 | nspec += 1 | 
|---|
| 444 |  | 
|---|
| 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, tabidx | 
|---|
| 454 |  | 
|---|
| 455 | return spec_array[0:nspec], tabidx | 
|---|
| 456 |  | 
|---|
| 457 |  | 
|---|
| 458 | @asaplog_post_dec | 
|---|
| 459 | def _setup_shift(self): | 
|---|
| 460 | ### define self.tables, self.signalShift, and self.imageShift | 
|---|
| 461 | 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." | 
|---|
| 467 |  | 
|---|
| 468 | byname = False | 
|---|
| 469 | #if not self.intables: | 
|---|
| 470 | if isinstance(self.intables[0], str): | 
|---|
| 471 | # A list of file name is given | 
|---|
| 472 | if not os.path.exists(self.intables[0]): | 
|---|
| 473 | asaplog.post() | 
|---|
| 474 | raise RuntimeError, "Could not find '%s'" % self.intables[0] | 
|---|
| 475 |  | 
|---|
| 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") | 
|---|
| 602 |  | 
|---|
| 603 | #     @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() | 
|---|