Changeset 2828


Ignore:
Timestamp:
07/30/13 15:18:19 (11 years ago)
Author:
Kana Sugimoto
Message:

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): sbseparator

Description: Removed obsolete codes and variables.


Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/python/sbseparator.py

    r2807 r2828  
    4040        self._separator = SBSeparator(infiles)
    4141
    42 ################## temp functions #################
    43 #     def rfft(self, invec):
    44 #         print "Python method cppRfft"
    45 #         self._separator._cpprfft(invec)
    46 #         return
    47 ################## temp functions #################
    4842
    4943    def set_frequency(self, baseif, freqtol, frame=""):
    … …  
    171165        self._separator.separate(outname)
    172166
    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)
    222        
    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."
    395        
    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]
    403            
    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)
    523 
    524 #     @asaplog_post_dec
    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()

  • trunk/src/STSideBandSep.cpp

    r2808 r2828  
    913913  os << "Threshold = " << threshold << ", Rejected channels = " << nreject << endl;
    914914};
    915 
    916 ////////////////////////////////////////////////////////////////////
    917 // void STSideBandSep::cpprfft(std::vector<float> invec)
    918 // {
    919 //   cout << "c++ method cpprfft" << endl;
    920 //   const unsigned int len = invec.size();
    921 //   Vector<Complex> carr(len/2+1, 0.);
    922 //   Vector<float> inarr = Vector<float>(invec);
    923 //   Vector <float> outarr(len, 0.);
    924 //   FFTServer<Float, Complex> fftsf, fftsi;
    925 //   fftsf.resize(IPosition(1, len), FFTEnums::REALTOCOMPLEX);
    926 //   fftsi.resize(IPosition(1, invec.size()), FFTEnums::COMPLEXTOREAL);
    927 //   cout << "Input float array (length = " << len << "):" << endl;
    928 //   for (uInt i = 0 ; i < len ; i++){
    929 //     cout << (i == 0 ? "( " : " ") << inarr[i] << (i == len-1 ? ")" : ",");
    930 //   }
    931 //   cout << endl;
    932 //   cout << "R->C transform" << endl;
    933 //   fftsf.fft0(carr, inarr, true);
    934 //   cout << "FFTed complex array (length = " << carr.size() << "):" << endl;
    935 //   for (uInt i = 0 ; i < carr.size() ; i++){
    936 //     cout << (i == 0 ? "( " : " ") << carr[i] << ( (i == carr.size()-1) ? ")" : "," );
    937 //   }
    938 //   cout << endl;
    939 //   cout << "C->R transform" << endl;
    940 //   fftsi.fft0(outarr, carr, false);
    941 //   cout << "invFFTed float array (length = " << outarr.size() << "):" << endl;
    942 //   for (uInt i = 0 ; i < outarr.size() ; i++){
    943 //     cout << (i == 0 ? "( " : " ") << outarr[i] << ( (i == outarr.size()-1) ? ")" : "," );
    944 //   }
    945 //   cout << endl;
    946 // };
    947 ////////////////////////////////////////////////////////////////////
    948915
    949916

  • trunk/src/STSideBandSep.h

    r2794 r2828  
    4242  virtual ~STSideBandSep();
    4343
    44   ///////////// temp functions //////////////////////
    45   //void cpprfft(std::vector<float> invec);
    46   //////////////////////////////////////////////////
    4744
    4845  /**
    … …  
    183180  Table::TableType tp_;
    184181  FFTServer<Float, Complex> fftsf, fftsi;
    185   // TEMPORAL member
    186   CountedPtr<Scantable> st0_;
    187182
    188183}; // class

  • trunk/src/python_STSideBandSep.cpp

    r2794 r2828  
    3636    .def( "separate", &STSideBandSep::separate )
    3737    //// temporal methods
    38     //.def( "_cpprfft", &STSideBandSep::cpprfft )
    3938    //.def( "solve_imgfreq", &STSideBandSep::solveImageFreqency )
    4039    //.def( "_get_asistb_from_scantb", &STSideBandSep::setScanTb0 )
Note: See TracChangeset for help on using the changeset viewer.