source: trunk/python/sbseparator.py @ 2821

Last change on this file since 2821 was 2807, checked in by Kana Sugimoto, 11 years ago

New Development: No (A bug fix for CASA 4.1)

JIRA Issue: Yes (CAS-4141/TRAC-290)

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): itself

Description: fixed a bug in sbseparator.set_lo1.


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