source: trunk/python/sbseparator.py@ 2700

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

New Development: No

JIRA Issue: Yes

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s):

Description: a bug fix and an update.
File size: 22.9 KB
Line 
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 """
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
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:
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
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', convsupport=1)
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")
514
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()
Note: See TracBrowser for help on using the repository browser.