source: trunk/python/sbseparator.py @ 2649

Last change on this file since 2649 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
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', 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.