source: trunk/python/sbseparator.py@ 2663

Last change on this file since 2663 was 2649, checked in by Kana Sugimoto, 12 years ago

New Development: No

JIRA Issue: Yes (CAS-4141/CSV-1526)

Ready for Test: Yes (* Fix for CASA 4.0 *)

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s):

Description:

For got to commit build and import code of sbseparator module.
Added known issues and warning messages.
File size: 23.0 KB
RevLine 
[2647]1import os, shutil
2import numpy
3import numpy.fft as FFT
4import math
5
6from asap.scantable import scantable
7from asap.parameters import rcParams
8from asap.logging import asaplog, asaplog_post_dec
9from asap.selector import selector
10from asap.asapgrid import asapgrid2
11#from asap._asap import sidebandsep
12
13class 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]))
350
351 asaplog.push("Regrid output scantable with cell = [%s, %s]" % \
352 (cellx, celly))
353 gridder.defineImage(nx=nx, ny=ny, cellx=cellx, celly=celly)
354 gridder.setFunc(func='box', width=1)
355 gridder.setWeight(weightType='uniform')
356 gridder.grid()
357 return gridder.getResult()
358
359
360 @asaplog_post_dec
361 def _get_specarray(self, polid=None, beamid=None, dir=None):
362 ntable = len(self.tables)
363 spec_array = numpy.zeros((ntable, self.nchan), numpy.float)
364 nspec = 0
365 asaplog.push("Start data selection by POL=%d, BEAM=%d, direction=[%f, %f]" % (polid, beamid, dir[0], dir[1]))
366 tabidx = []
367 for itab in range(ntable):
368 tab = self.tables[itab]
369 # Select rows by POLNO and BEAMNO
370 try:
371 tab.set_selection(pols=[polid], beams=[beamid])
372 if tab.nrow() > 0: tabidx.append(itab)
373 except: # no selection
374 asaplog.post()
375 asaplog.push("table %d - No spectrum ....skipping the table" % (itab))
376 asaplog.post("WARN")
377 continue
378
379 # Select rows by direction
380 spec = numpy.zeros(self.nchan, numpy.float)
381 selrow = []
382 for irow in range(tab.nrow()):
383 currdir = tab.get_directionval(irow)
384 if (abs(currdir[0] - dir[0]) > self.dirtol[0]) or \
385 (abs(currdir[1] - dir[1]) > self.dirtol[1]):
386 continue
387 selrow.append(irow)
388 if len(selrow) == 0:
389 asaplog.post()
390 asaplog.push("table %d - No spectrum ....skipping the table" % (itab))
391 asaplog.post("WARN")
392 continue
393 else:
394 seltab = tab.copy()
395 seltab.set_selection(selector(rows=selrow))
396
397 if tab.nrow() > 1:
398 asaplog.push("table %d - More than a spectrum selected. averaging rows..." % (itab))
399 tab = seltab.average_time(scanav=False, weight="tintsys")
400 else:
401 tab = seltab
402
403 spec_array[nspec] = tab._getspectrum()
404 nspec += 1
405
406 if nspec != ntable:
407 asaplog.post()
408 #asaplog.push("Some tables has no spectrum with POL=%d BEAM=%d. averaging rows..." % (polid, beamid))
409 asaplog.push("Could not find corresponding rows in some tables.")
410 asaplog.push("Number of spectra selected = %d (table: %d)" % (nspec, ntable))
411 if nspec < 2:
412 asaplog.push("At least 2 spectra are necessary for convolution")
413 asaplog.post("ERROR")
414 return False, tabidx
415
416 return spec_array[0:nspec], tabidx
417
418
419 @asaplog_post_dec
420 def _setup_shift(self):
421 ### define self.tables, self.signalShift, and self.imageShift
422 if not self.intables:
423 asaplog.post()
424 raise RuntimeError, "Input data is not defined."
425 #if self.baseif < 0:
426 # asaplog.post()
427 # raise RuntimeError, "Reference IFNO is not defined."
428
429 byname = False
430 #if not self.intables:
431 if isinstance(self.intables[0], str):
432 # A list of file name is given
433 if not os.path.exists(self.intables[0]):
434 asaplog.post()
435 raise RuntimeError, "Could not find '%s'" % self.intables[0]
436
437 stab = scantable(self.intables[0],average=False)
438 ntab = len(self.intables)
439 byname = True
440 else:
441 stab = self.intables[0]
442 ntab = len(self.intables)
443
444 if len(stab.getbeamnos()) > 1:
445 asaplog.post()
446 asaplog.push("Mult-beam data is not supported by this module.")
447 asaplog.post("ERROR")
448 return
449
450 valid_ifs = stab.getifnos()
451 if self.baseif < 0:
452 self.baseif = valid_ifs[0]
453 asaplog.post()
454 asaplog.push("IFNO is not selected. Using the first IF in the first scantable. Reference IFNO = %d" % (self.baseif))
455
456 if not (self.baseif in valid_ifs):
457 asaplog.post()
458 errmsg = "IF%d does not exist in the first scantable" % \
459 self.baseif
460 raise RuntimeError, errmsg
461
462 asaplog.push("Start selecting tables and IFNOs to solve.")
463 asaplog.push("Cheching frequency of the reference IF")
464 unit_org = stab.get_unit()
465 coord = stab._getcoordinfo()
466 frame_org = coord[1]
467 stab.set_unit("Hz")
468 if len(self.freqframe) > 0:
469 stab.set_freqframe(self.freqframe)
470 stab.set_selection(ifs=[self.baseif])
471 spx = stab._getabcissa()
472 stab.set_selection()
473 basech0 = spx[0]
474 baseinc = spx[1]-spx[0]
475 self.nchan = len(spx)
476 if isinstance(self.freqtol, float):
477 vftol = abs(baseinc * self.freqtol)
478 self.freqtol = dict(value=vftol, unit="Hz")
479 else:
480 vftol = abs(self.freqtol['value'])
481 inctol = abs(baseinc/float(self.nchan))
482 asaplog.push("Reference frequency setup (Table = 0, IFNO = %d): nchan = %d, chan0 = %f Hz, incr = %f Hz" % (self.baseif, self.nchan, basech0, baseinc))
483 asaplog.push("Allowed frequency tolerance = %f Hz ( %f channels)" % (vftol, vftol/baseinc))
484 poltype0 = stab.poltype()
485
486 self.tables = []
487 self.signalShift = []
488 if self.dsbmode:
489 self.imageShift = []
490
491 for itab in range(ntab):
492 asaplog.push("Table %d:" % itab)
493 tab_selected = False
494 if itab > 0:
495 if byname:
496 stab = scantable(self.intables[itab],average=False)
497 self.intables.append(stab)
498 else:
499 stab = self.intables[itab]
500 unit_org = stab.get_unit()
501 coord = stab._getcoordinfo()
502 frame_org = coord[1]
503 stab.set_unit("Hz")
504 if len(self.freqframe) > 0:
505 stab.set_freqframe(self.freqframe)
506
507 # Check POLTYPE should be equal to the first table.
508 if stab.poltype() != poltype0:
509 asaplog.post()
510 raise Exception, "POLTYPE should be equal to the first table."
511 # Multiple beam data may not handled properly
512 if len(stab.getbeamnos()) > 1:
513 asaplog.post()
514 asaplog.push("table contains multiple beams. It may not be handled properly.")
515 asaplog.push("WARN")
516
517 for ifno in stab.getifnos():
518 stab.set_selection(ifs=[ifno])
519 spx = stab._getabcissa()
520 if (len(spx) != self.nchan) or \
521 (abs(spx[0]-basech0) > vftol) or \
522 (abs(spx[1]-spx[0]-baseinc) > inctol):
523 continue
524 tab_selected = True
525 seltab = stab.copy()
526 seltab.set_unit(unit_org)
527 seltab.set_freqframe(frame_org)
528 self.tables.append(seltab)
529 self.signalShift.append((spx[0]-basech0)/baseinc)
530 if self.dsbmode:
531 self.imageShift.append(-self.signalShift[-1])
532 asaplog.push("- IF%d selected: sideband shift = %16.12e channels" % (ifno, self.signalShift[-1]))
533 stab.set_selection()
534 stab.set_unit(unit_org)
535 stab.set_freqframe(frame_org)
536 if not tab_selected:
537 asaplog.post()
538 asaplog.push("No data selected in Table %d" % itab)
539 asaplog.post("WARN")
540
541 asaplog.push("Total number of IFs selected = %d" % len(self.tables))
542 if len(self.tables) < 2:
543 asaplog.post()
544 raise RuntimeError, "At least 2 IFs are necessary for convolution!"
545
546 if not self.dsbmode and len(self.imageShift) != len(self.signalShift):
547 asaplog.post()
548 errmsg = "User defined channel shift of image sideband has %d elements, while selected IFNOs are %d" % (len(self.imageShift), len(self.signalShift))
549 raise RuntimeError, errmsg
550
551 self.signalShift = numpy.array(self.signalShift)
552 self.imageShift = numpy.array(self.imageShift)
553 self.nshift = len(self.tables)
554
555 @asaplog_post_dec
556 def _preprocess_tables(self):
557 ### temporary method to preprocess data
558 ### Do time averaging for now.
559 for itab in range(len(self.tables)):
560 self.tables[itab] = self.tables[itab].average_time(scanav=False, weight="tintsys")
561
562
563# def save(self, outfile, outform="ASAP", overwrite=False):
564# if not overwrite and os.path.exists(outfile):
565# raise RuntimeError, "Output file '%s' already exists" % outfile
566#
567# #self.separator._save(outfile, outform)
568
569# def done(self):
570# self.close()
571
572# def close(self):
573# pass
574# #del self.separator
575
576
577
578########################################################################
579 def _Deconvolution(self, data_array, shift, threshold=0.00000001):
580 FObs = []
581 Reject = 0
582 nshift, nchan = data_array.shape
583 nspec = nshift*(nshift-1)/2
584 ifftObs = numpy.zeros((nspec, nchan), numpy.float)
585 for i in range(nshift):
586 F = FFT.fft(data_array[i])
587 FObs.append(F)
588 z = 0
589 for i in range(nshift):
590 for j in range(i+1, nshift):
591 Fobs = (FObs[i]+FObs[j])/2.0
592 dX = (shift[j]-shift[i])*2.0*math.pi/float(self.nchan)
593 #print 'dX,i,j=',dX,i,j
594 for k in range(1,self.nchan):
595 if math.fabs(math.sin(dX*k)) > threshold:
596 Fobs[k] += ((FObs[i][k]-FObs[j][k])/2.0/(1.0-math.cos(dX*k))*math.sin(dX*k))*1.0j
597 else: Reject += 1
598 ifftObs[z] = FFT.ifft(Fobs)
599 z += 1
600 print 'Threshold=%s Reject=%d' % (threshold, Reject)
601 return ifftObs
602
603 def _combineResult(self, ifftObs):
604 nspec = len(ifftObs)
605 sum = ifftObs[0]
606 for i in range(1,nspec):
607 sum += ifftObs[i]
608 return(sum/float(nspec))
609
610 def _subtractOtherSide(self, data_array, shift, Data):
611 sum = numpy.zeros(len(Data), numpy.float)
612 numSP = len(data_array)
613 for i in range(numSP):
614 SPsub = data_array[i] - Data
615 sum += self._shiftSpectrum(SPsub, -shift[i])
616 return(sum/float(numSP))
617
618 def _shiftSpectrum(self, data, Shift):
619 Out = numpy.zeros(self.nchan, numpy.float)
620 w2 = Shift % 1
621 w1 = 1.0 - w2
622 for i in range(self.nchan):
623 c1 = int((Shift + i) % self.nchan)
624 c2 = (c1 + 1) % self.nchan
625 Out[c1] += data[i] * w1
626 Out[c2] += data[i] * w2
627 return Out.copy()
Note: See TracBrowser for help on using the repository browser.