Changeset 2238


Ignore:
Timestamp:
07/21/11 12:17:06 (13 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

Improved performance.
Speed up getDirection.


Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/MSFiller.cpp

    r2237 r2238  
    263263    pointtab = MSPointing( pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort("TIME") ) ;
    264264  }
     265  ROScalarColumn<Double> ptcol( pointtab, "TIME" ) ;
     266  ROArrayColumn<Double> pdcol( pointtab, "DIRECTION" ) ;
    265267  String stationName = asString( "STATION", antenna_, anttab, tpoolr ) ;
    266268  String antennaName = asString( "NAME", antenna_, anttab, tpoolr ) ;
    … …  
    651653                String refString ;
    652654                if ( getPt_ )
    653                   diridx = getDirection( diridx, dir, azel, scanrate, pointtab, mTimeB[irow], mp ) ;
     655                  //diridx = getDirection( diridx, dir, azel, scanrate, pointtab, mTimeB[irow], mp ) ;
     656                  diridx = getDirection( diridx, dir, azel, scanrate, ptcol, pdcol, mTimeB[irow], mp ) ;
    654657                else
    655658                  getSourceDirection( dir, azel, scanrate, mTimeB[irow], mp, delayDir ) ;
    … …  
    13991402}
    14001403
    1401 uInt MSFiller::getDirection( uInt idx, Vector<Double> &dir, Vector<Double> &srate, String &ref, MSPointing &tab, Double t )
     1404uInt MSFiller::getDirection( uInt idx,
     1405                             Vector<Double> &dir,
     1406                             Vector<Double> &srate,
     1407                             String &ref,
     1408                             ROScalarColumn<Double> &tcol,
     1409                             ROArrayColumn<Double> &dcol,
     1410                             Double t )
    14021411{
    14031412  //double startSec = gettimeofday_sec() ;
    14041413  //os_ << "start MSFiller::getDirection1() startSec=" << startSec << LogIO::POST ;
     1414  //double time0 = gettimeofday_sec() ;
     1415  //os_ << "start getDirection 1st stage startSec=" << time0 << LogIO::POST ;
    14051416  // assume that cols is sorted by TIME
    14061417  Bool doInterp = False ;
    1407   //uInt nrow = cols.nrow() ;
    1408   uInt nrow = tab.nrow() ;
     1418  uInt nrow = tcol.nrow() ;
    14091419  if ( nrow == 0 )
    14101420    return 0 ;
    1411   ROScalarMeasColumn<MEpoch> tcol( tab, "TIME" ) ;
    1412   ROArrayMeasColumn<MDirection> dmcol( tab, "DIRECTION" ) ;
    1413   ROArrayColumn<Double> dcol( tab, "DIRECTION" ) ;
     1421  TableRecord trec = tcol.keywordSet() ;
     1422  String tUnit = trec.asArrayString( "QuantumUnits" ).data()[0] ;
     1423  Double factor = 1.0 ;
     1424  if ( tUnit == "d" )
     1425    factor = 86400.0 ;
    14141426  // binary search if idx == 0
    14151427  if ( idx == 0 ) {
    1416     uInt nrowb = 75000 ;
     1428    uInt nrowb = 1000 ;
    14171429    if ( nrow > nrowb ) {
    14181430      uInt nblock = nrow / nrowb + 1 ;
    … …  
    14201432        uInt high = min( nrowb, nrow-iblock*nrowb ) ;
    14211433
    1422         if ( tcol( high-1 ).get( "s" ).getValue() < t ) {
     1434        if ( tcol( high-1 ) * factor < t ) {
    14231435          idx = iblock * nrowb ;
    14241436          continue ;
    14251437        }
    14261438
    1427         Vector<MEpoch> tarr( high ) ;
    1428         for ( uInt irow = 0 ; irow < high ; irow++ ) {
    1429           tarr[irow] = tcol( iblock*nrowb+irow ) ;
    1430         }
     1439        RefRows refrows( iblock*nrowb, iblock*nrowb+high, 1 ) ;
     1440        Vector<Double> tarr = tcol.getColumnCells( refrows ) ;
     1441        if ( tUnit == "d" )
     1442          tarr *= factor ;
    14311443
    14321444        uInt bidx = binarySearch( tarr, t ) ;
    … …  
    14371449    }
    14381450    else {
    1439       Vector<MEpoch> tarr( nrow ) ;
    1440       for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
    1441         tarr[irow] = tcol( irow ) ;
    1442       }
     1451      Vector<Double> tarr = tcol.getColumn() ;
     1452      if ( tUnit == "d" )
     1453        tarr *= factor ;
    14431454      idx = binarySearch( tarr, t ) ;
    14441455    }
    14451456  }
     1457  //double time1 = gettimeofday_sec() ;
     1458  //os_ << "end getDirection 1st stage endSec=" << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    14461459  // ensure that tcol(idx) < t
    14471460  //os_ << "tcol(idx) = " << tcol(idx).get("s").getValue() << " t = " << t << " diff = " << tcol(idx).get("s").getValue()-t << endl ;
    1448   while ( tcol(idx).get("s").getValue() > t && idx > 0 )
     1461  //time0 = gettimeofday_sec() ;
     1462  //os_ << "start getDirection 2nd stage startSec=" << time0 << LogIO::POST ;
     1463  while( tcol( idx ) * factor > t && idx > 0 )
    14491464    idx-- ;
    14501465  //os_ << "idx = " << idx << LogIO::POST ;
    … …  
    14521467  // index search
    14531468  for ( uInt i = idx ; i < nrow ; i++ ) {
    1454     Double tref = tcol( i ).get( "s" ).getValue() ;
     1469    //Double tref = tcol( i ).get( "s" ).getValue() ;
     1470    Double tref = tcol( i ) * factor ;
    14551471    if ( tref == t ) {
    14561472      idx = i ;
    … …  
    14711487    }
    14721488  }
     1489  //time1 = gettimeofday_sec() ;
     1490  //os_ << "end getDirection 2nd stage endSec=" << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    14731491  //os_ << "searched idx = " << idx << LogIO::POST ;
    14741492
     1493  //time0 = gettimeofday_sec() ;
     1494  //os_ << "start getDirection 3rd stage startSec=" << time0 << LogIO::POST ;
    14751495  //os_ << "dmcol(idx).shape() = " << dmcol(idx).shape() << LogIO::POST ;
    1476   IPosition ip( dmcol(idx).shape().nelements(), 0 ) ;
     1496  //IPosition ip( dmcol(idx).shape().nelements(), 0 ) ;
     1497  IPosition ip( dcol(idx).shape().nelements(), 0 ) ;
    14771498  //os_ << "ip = " << ip << LogIO::POST ;
    1478   ref = dmcol(idx)(ip).getRefString() ;
     1499  //ref = dmcol(idx)(ip).getRefString() ;
     1500  trec = dcol.keywordSet() ;
     1501  Record rec = trec.asRecord( "MEASINFO" ) ;
     1502  ref = rec.asString( "Ref" ) ;
    14791503  //os_ << "ref = " << ref << LogIO::POST ;
    14801504  if ( doInterp ) {
    14811505    //os_ << "do interpolation" << LogIO::POST ;
    14821506    //os_ << "dcol(idx).shape() = " << dcol(idx).shape() << LogIO::POST ;
    1483     Double tref0 = tcol(idx).get("s").getValue() ;
    1484     Double tref1 = tcol(idx+1).get("s").getValue() ;
     1507//     Double tref0 = tcol(idx).get("s").getValue() ;
     1508//     Double tref1 = tcol(idx+1).get("s").getValue() ;
     1509    Double tref0 = tcol(idx) * factor ;
     1510    Double tref1 = tcol(idx+1) * factor ;
    14851511    Matrix<Double> mdir0 = dcol( idx ) ;
    14861512    Matrix<Double> mdir1 = dcol( idx+1 ) ;
    … …  
    15081534  }
    15091535
     1536  //time1 = gettimeofday_sec() ;
     1537  //os_ << "end getDirection 3rd stage endSec=" << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    15101538  //double endSec = gettimeofday_sec() ;
    15111539  //os_ << "end MSFiller::getDirection1() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    … …  
    15331561    else if ( t > target )
    15341562      high = idx - 1 ;
    1535     else
     1563    else {
    15361564      return idx ;
     1565    }
    15371566  }
    15381567
    … …  
    15401569
    15411570  return idx ;
    1542  
     1571}
     1572
     1573uInt MSFiller::binarySearch( Vector<Double> &timeList, Double target )
     1574{
     1575  Int low = 0 ;
     1576  Int high = timeList.nelements() ;
     1577  uInt idx = 0 ;
     1578
     1579  while ( low <= high ) {
     1580    idx = (Int)( 0.5 * ( low + high ) ) ;
     1581    Double t = timeList[idx] ;
     1582    if ( t < target )
     1583      low = idx + 1 ;
     1584    else if ( t > target )
     1585      high = idx - 1 ;
     1586    else {
     1587      return idx ;
     1588    }
     1589  }
     1590
     1591  idx = max( 0, min( low, high ) ) ;
     1592
     1593  return idx ;
    15431594}
    15441595
    … …  
    15841635    ROArrayColumn<Bool> mFlagCol( tab, "FLAG" ) ;
    15851636    ROArrayColumn<Complex> mDataCol( tab, "DATA" ) ;
    1586     for ( Int irow = 0 ; irow < nrow ; irow++ ) {
    1587       Bool crossOK = False ;
    1588       Matrix<Complex> mSp = mDataCol( irow ) ;
    1589       Matrix<Bool> mFl = mFlagCol( irow ) ;
    1590       Matrix<Float> spxy = sp.xyPlane( irow ) ;
    1591       Matrix<Bool> flxy = fl.xyPlane( irow ) ;
    1592       for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
    1593         if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::YX
    1594              || corrtype[ipol] == Stokes::RL || corrtype[ipol] == Stokes::LR ) {
    1595           if ( !crossOK ) {
     1637    if ( npol < 3 ) {
     1638//       Cube<Complex> tmp = mDataCol.getColumn() ;
     1639//       Float *sp_p = sp.data() ;
     1640//       Complex *tmp_p = tmp.data() ;
     1641//       for ( uInt i = 0 ; i < sp.nelements() ; i++ )
     1642//         sp_p[i] = tmp_p[i].real() ;
     1643      Cube<Float> tmp = ComplexToReal( mDataCol.getColumn() ) ;
     1644      IPosition start( 3, 0, 0, 0 ) ;
     1645      IPosition end( 3, 2*npol-1, nchan-1, nrow-1 ) ;
     1646      IPosition inc( 3, 2, 1, 1 ) ;
     1647      sp = tmp( start, end, inc ) ;
     1648//       sp = tmp( Slice( 0, npol, 2 ),
     1649//                 Slice( 0, nchan, 1 ),
     1650//                 Slice( 0, nrow, 1 ) ) ;
     1651      //sp = real( mDataCol.getColumn() ) ;
     1652      fl = mFlagCol.getColumn() ;
     1653    }
     1654    else {
     1655      for ( Int irow = 0 ; irow < nrow ; irow++ ) {
     1656        Bool crossOK = False ;
     1657        Matrix<Complex> mSp = mDataCol( irow ) ;
     1658        Matrix<Bool> mFl = mFlagCol( irow ) ;
     1659        Matrix<Float> spxy = sp.xyPlane( irow ) ;
     1660        Matrix<Bool> flxy = fl.xyPlane( irow ) ;
     1661        for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
     1662          if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::YX
     1663               || corrtype[ipol] == Stokes::RL || corrtype[ipol] == Stokes::LR ) {
     1664            if ( !crossOK ) {
     1665              spxy.row( ipol ) = real( mSp.row( ipol ) ) ;
     1666              flxy.row( ipol ) = mFl.row( ipol ) ;
     1667              if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::RL ) {
     1668                spxy.row( ipol+1 ) = imag( mSp.row( ipol ) ) ;
     1669                flxy.row( ipol+1 ) = mFl.row( ipol ) ;
     1670              }                       
     1671              else {
     1672                spxy.row( ipol+1 ) = imag( conj( mSp.row( ipol ) ) ) ;
     1673                flxy.row( ipol+1 ) = mFl.row( ipol ) ;
     1674              }
     1675              crossOK = True ;
     1676            }
     1677          }
     1678          else {
    15961679            spxy.row( ipol ) = real( mSp.row( ipol ) ) ;
    15971680            flxy.row( ipol ) = mFl.row( ipol ) ;
    1598             if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::RL ) {
    1599               spxy.row( ipol+1 ) = imag( mSp.row( ipol ) ) ;
    1600               flxy.row( ipol+1 ) = mFl.row( ipol ) ;
    1601             }                       
    1602             else {
    1603               spxy.row( ipol+1 ) = imag( conj( mSp.row( ipol ) ) ) ;
    1604               flxy.row( ipol+1 ) = mFl.row( ipol ) ;
    1605             }
    1606             crossOK = True ;
    16071681          }
    1608         }
    1609         else {
    1610           spxy.row( ipol ) = real( mSp.row( ipol ) ) ;
    1611           flxy.row( ipol ) = mFl.row( ipol ) ;
    16121682        }
    16131683      }
    … …  
    16221692                             Vector<Double> &azel,
    16231693                             Vector<Double> &srate,
    1624                              MSPointing &ptab,
     1694                             ROScalarColumn<Double> &ptcol,
     1695                             ROArrayColumn<Double> &pdcol,
    16251696                             MEpoch &t,
    16261697                             MPosition &antpos )
    … …  
    16301701  String refString ;
    16311702  MDirection::Types dirType ;
    1632   uInt diridx = getDirection( idx, dir, srate, refString, ptab, t.get("s").getValue() ) ;
     1703  //uInt diridx = getDirection( idx, dir, srate, refString, ptab, t.get("s").getValue() ) ;
     1704  uInt diridx = getDirection( idx, dir, srate, refString, ptcol, pdcol, t.get("s").getValue() ) ;
    16331705  MDirection::getType( dirType, refString ) ;
    16341706  MeasFrame mf( t, antpos ) ;

  • trunk/src/MSFiller.h

    r2237 r2238  
    9090
    9191  // get direction for DIRECTION, AZIMUTH, and ELEVATION columns
    92   casa::uInt getDirection( casa::uInt idx, casa::Vector<casa::Double> &dir, casa::Vector<casa::Double> &srate, casa::String &ref, casa::MSPointing &tab, casa::Double t ) ;
    93   casa::uInt getDirection( casa::uInt idx, casa::Vector<casa::Double> &dir, casa::Vector<casa::Double> &azel, casa::Vector<casa::Double> &srate, casa::MSPointing &tab, casa::MEpoch &t, casa::MPosition &antpos ) ;
    94   void getSourceDirection( casa::Vector<casa::Double> &dir, casa::Vector<casa::Double> &azel, casa::Vector<casa::Double> &srate, casa::MEpoch &t, casa::MPosition &antpos, casa::Vector<casa::MDirection> &srcdir ) ;
     92  casa::uInt getDirection( casa::uInt idx,
     93                           casa::Vector<casa::Double> &dir,
     94                           casa::Vector<casa::Double> &srate,
     95                           casa::String &ref,
     96                           casa::ROScalarColumn<casa::Double> &ptcol,
     97                           casa::ROArrayColumn<casa::Double> &pdcol,
     98                           casa::Double t ) ;
     99  casa::uInt getDirection( casa::uInt idx,
     100                           casa::Vector<casa::Double> &dir,
     101                           casa::Vector<casa::Double> &azel,
     102                           casa::Vector<casa::Double> &srate,
     103                           casa::ROScalarColumn<casa::Double> &ptcol,
     104                           casa::ROArrayColumn<casa::Double> &pdcol,
     105                           casa::MEpoch &t, casa::MPosition &antpos ) ;
     106  void getSourceDirection( casa::Vector<casa::Double> &dir,
     107                           casa::Vector<casa::Double> &azel,
     108                           casa::Vector<casa::Double> &srate,
     109                           casa::MEpoch &t,
     110                           casa::MPosition &antpos,
     111                           casa::Vector<casa::MDirection> &srcdir ) ;
    95112
    96113  // create key for TCAL table
    … …  
    99116  // binary search
    100117  casa::uInt binarySearch( casa::Vector<casa::MEpoch> &timeList, casa::Double target ) ;
     118  casa::uInt binarySearch( casa::Vector<casa::Double> &timeList, casa::Double target ) ;
    101119
    102120  // tool for HPC
Note: See TracChangeset for help on using the changeset viewer.