Changeset 2794 for trunk/python


Ignore:
Timestamp:
03/15/13 18:23:04 (11 years ago)
Author:
Kana Sugimoto
Message:

New Development: Yes

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): sbseparator module

Description: Pushed python codes down to cpp.


File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/python/sbseparator.py

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