Changeset 1677 for branches/alma


Ignore:
Timestamp:
01/27/10 18:29:48 (14 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: Yes CAS-1823

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: The mode parameter is added to scantable.scale() method.

Test Programs: s = sd.scantable('yourfile',False)

factor = []
for i in range(s.nrow()):

factor.append(i)

s2 = s + factor

Put in Release Notes: Yes

Module(s): -

Description: Describe your changes here...

Basic operations (addition, subtraction, multiplication, division)
of scantable with one dimensional list are implemented.
Size of list operand should be equal to either number of spectral channel
or number of row. In the former case, the list is operated as
channel-by-channel manner, while it is operated as row-by-row manner
in the latter case.
If number of spectral channel is equal to number of row, row-by-row
operation will be done.

The user is able to select operation mode (channel-by-channel or row-by-row)
manually by using lower level function, stmath.arrayop().

The scantable.scale() method is updated to allow list scaling factor.
Scaling is done in channel-by-channel manner if mode is set to 'channel',
while in row-by-row manner if mode is set to 'row'.


Location:
branches/alma
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • branches/alma/python/scantable.py

    r1676 r1677  
    18681868            return s
    18691869
    1870     def scale(self, factor, tsys=True, insitu=None):
     1870    def scale(self, factor, tsys=True, mode='channel', insitu=None):
    18711871        """
    18721872        Return a scan where all spectra are scaled by the give 'factor'
    18731873        Parameters:
    1874             factor:      the scaling factor
     1874            factor:      the scaling factor (float or 1D float list)
    18751875            insitu:      if False a new scantable is returned.
    18761876                         Otherwise, the scaling is done in-situ
    … …  
    18781878            tsys:        if True (default) then apply the operation to Tsys
    18791879                         as well as the data
     1880            mode:        operation mode for list scaling factor. Possible
     1881                         values are 'channel' (channel-by-channel) or
     1882                         'row' (row-by-row).
    18801883        """
    18811884        if insitu is None: insitu = rcParams['insitu']
    18821885        self._math._setinsitu(insitu)
    18831886        varlist = vars()
    1884         s = scantable(self._math._unaryop(self, factor, "MUL", tsys))
     1887        import numpy
     1888        if type(factor) == list:
     1889            if mode == 'row':
     1890                s = scantable( self._math._arrayop( self.copy(), factor, "MUL", tsys, 'row' ) )
     1891            elif mode == 'channel':
     1892                s = scantable( self._math._arrayop( self.copy(), factor, "MUL", tsys, 'channel' ) )
     1893            else:
     1894                asaplog.push( 'scantable.scale(): Unknown operation mode. No scaling.' )
     1895                print_log()
     1896                s = self.copy()
     1897        else:
     1898            s = scantable(self._math._unaryop(self.copy(), factor, "MUL", tsys))
    18851899        s._add_history("scale", varlist)
    18861900        print_log()
    … …  
    20272041        s = None
    20282042        if isinstance(other, scantable):
    2029             s = scantable(self._math._binaryop(self, other, "ADD"))
     2043            s = scantable(self._math._binaryop(self.copy(), other, "ADD"))
    20302044        elif isinstance(other, float):
    2031             s = scantable(self._math._unaryop(self, other, "ADD", False))
    2032         else:
    2033             raise TypeError("Other input is not a scantable or float value")
     2045            s = scantable(self._math._unaryop(self.copy(), other, "ADD", False))
     2046        elif isinstance(other, list):
     2047            if len(other) == self.nrow():
     2048                s = scantable(self._math._arrayop(self.copy(), other, "ADD", False, 'row'))
     2049            elif len(other) == self.nchan(0):
     2050                s = scantable(self._math._arrayop(self.copy(), other, "ADD", False, 'channel'))
     2051            else:
     2052                casalog.post( 'Length of list operand should equal to either scantable.nchan() or scantable.nrow().', 'ERROR' )
     2053        else:
     2054            raise TypeError("Other input is not a scantable or float value or float list")
    20342055        s._add_history("operator +", varlist)
    20352056        print_log()
    … …  
    20432064        s = None
    20442065        if isinstance(other, scantable):
    2045             s = scantable(self._math._binaryop(self, other, "SUB"))
     2066            s = scantable(self._math._binaryop(self.copy(), other, "SUB"))
    20462067        elif isinstance(other, float):
    2047             s = scantable(self._math._unaryop(self, other, "SUB", False))
    2048         else:
    2049             raise TypeError("Other input is not a scantable or float value")
     2068            s = scantable(self._math._unaryop(self.copy(), other, "SUB", False))
     2069        elif isinstance(other, list):
     2070            if len(other) == self.nrow():
     2071                s = scantable(self._math._arrayop(self.copy(), other, "SUB", False, 'row'))
     2072            elif len(other) == self.nchan(0):
     2073                s = scantable(self._math._arrayop(self.copy(), other, "SUB", False, 'channel'))
     2074            else:
     2075                casalog.post( 'Length of list operand should equal to either scantable.nchan() or scantable.nrow().', 'ERROR' )
     2076        else:
     2077            raise TypeError("Other input is not a scantable or float value or float list")
    20502078        s._add_history("operator -", varlist)
    20512079        print_log()
    … …  
    20592087        s = None
    20602088        if isinstance(other, scantable):
    2061             s = scantable(self._math._binaryop(self, other, "MUL"))
     2089            s = scantable(self._math._binaryop(self.copy(), other, "MUL"))
    20622090        elif isinstance(other, float):
    2063             s = scantable(self._math._unaryop(self, other, "MUL", False))
    2064         else:
    2065             raise TypeError("Other input is not a scantable or float value")
     2091            s = scantable(self._math._unaryop(self.copy(), other, "MUL", False))
     2092        elif isinstance(other, list):
     2093            if len(other) == self.nrow():
     2094                s = scantable(self._math._arrayop(self.copy(), other, "MUL", False, 'row'))
     2095            elif len(other) == self.nchan(0):
     2096                s = scantable(self._math._arrayop(self.copy(), other, "MUL", False, 'channel'))
     2097            else:
     2098                casalog.post( 'Length of list operand should equal to either scantable.nchan() or scantable.nrow().', 'ERROR' )
     2099        else:
     2100            raise TypeError("Other input is not a scantable or float value or float list")
    20662101        s._add_history("operator *", varlist)
    20672102        print_log()
    20682103        return s
    2069 
    20702104
    20712105    def __div__(self, other):
    … …  
    20762110        s = None
    20772111        if isinstance(other, scantable):
    2078             s = scantable(self._math._binaryop(self, other, "DIV"))
     2112            s = scantable(self._math._binaryop(self.copy(), other, "DIV"))
    20792113        elif isinstance(other, float):
    20802114            if other == 0.0:
    20812115                raise ZeroDivisionError("Dividing by zero is not recommended")
    2082             s = scantable(self._math._unaryop(self, other, "DIV", False))
    2083         else:
    2084             raise TypeError("Other input is not a scantable or float value")
     2116            s = scantable(self._math._unaryop(self.copy(), other, "DIV", False))
     2117        elif isinstance(other, list):
     2118            if len(other) == self.nrow():
     2119                s = scantable(self._math._array2dop(self.copy(), other, "DIV", False, 'row'))
     2120            elif len(other) == self.nchan(0):
     2121                s = scantable(self._math._arrayop(self.copy(), other, "DIV", False, 'channel'))
     2122            else:
     2123                casalog.post( 'Length of list operand should equal to either scantable.nchan() or scantable.nrow().', 'ERROR' )
     2124        else:
     2125            raise TypeError("Other input is not a scantable or float value or float list")
    20852126        s._add_history("operator /", varlist)
    20862127        print_log()

  • branches/alma/src/STMath.cpp

    r1675 r1677  
    497497      if ( tsys ) {
    498498        ts += val;
     499        tsysCol.put(i, ts);
     500      }
     501    }
     502  }
     503  return out;
     504}
     505
     506CountedPtr< Scantable > STMath::arrayOperate( const CountedPtr< Scantable >& in,
     507                                              const std::vector<float> val,
     508                                              const std::string& mode,
     509                                              bool tsys,
     510                                              const std::string& op )
     511{
     512  CountedPtr< Scantable > out ;
     513  if ( op == "channel" ) {
     514    out = arrayOperateChannel( in, val, mode, tsys ) ;
     515  }
     516  else if ( op == "row" ) {
     517    out = arrayOperateRow( in, val, mode, tsys ) ;
     518  }
     519  else {
     520    throw( AipsError( "Unknown array operation mode." ) ) ;
     521  }
     522  return out ;
     523}
     524
     525CountedPtr< Scantable > STMath::arrayOperateChannel( const CountedPtr< Scantable >& in,
     526                                                     const std::vector<float> val,
     527                                                     const std::string& mode,
     528                                                     bool tsys )
     529{
     530  // conformity of SPECTRA and TSYS
     531  if ( tsys ) {
     532    TableIterator titer(in->table(), "IFNO");
     533    while ( !titer.pastEnd() ) {
     534      ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
     535      ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
     536      Array<Float> spec = specCol.getColumn() ;
     537      Array<Float> ts = tsysCol.getColumn() ;
     538      if ( !spec.conform( ts ) ) {
     539        throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
     540      }
     541      titer.next() ;
     542    }
     543  }
     544
     545  // check if all spectra in the scantable have the same number of channel
     546  vector<uInt> nchans;
     547  vector<uInt> ifnos = in->getIFNos() ;
     548  for ( uInt i = 0 ; i < ifnos.size() ; i++ ) {
     549    nchans.push_back( in->nchan( ifnos[i] ) ) ;
     550  }
     551  Vector<uInt> mchans( nchans ) ;
     552  if ( anyNE( mchans, mchans[0] ) ) {
     553    throw( AipsError("All spectra in the input scantable must have the same number of channel for vector operation." ) ) ;
     554  }
     555
     556  // check if vector size is equal to nchan
     557  Vector<Float> fact( val ) ;
     558  if ( fact.nelements() != mchans[0] ) {
     559    throw( AipsError("Vector size must be same as number of channel.") ) ;
     560  }
     561
     562  // check divided by zero
     563  if ( ( mode == "DIV" ) && anyEQ( fact, (float)0.0 ) ) {
     564    throw( AipsError("Divided by zero is not recommended." ) ) ;
     565  }
     566
     567  CountedPtr< Scantable > out = getScantable(in, false);
     568  Table& tab = out->table();
     569  ArrayColumn<Float> specCol(tab,"SPECTRA");
     570  ArrayColumn<Float> tsysCol(tab,"TSYS");
     571  for (uInt i=0; i<tab.nrow(); ++i) {
     572    Vector<Float> spec;
     573    Vector<Float> ts;
     574    specCol.get(i, spec);
     575    tsysCol.get(i, ts);
     576    if (mode == "MUL" || mode == "DIV") {
     577      if (mode == "DIV") fact = (float)1.0 / fact;
     578      spec *= fact;
     579      specCol.put(i, spec);
     580      if ( tsys ) {
     581        ts *= fact;
     582        tsysCol.put(i, ts);
     583      }
     584    } else if ( mode == "ADD"  || mode == "SUB") {
     585      if (mode == "SUB") fact *= (float)-1.0 ;
     586      spec += fact;
     587      specCol.put(i, spec);
     588      if ( tsys ) {
     589        ts += fact;
     590        tsysCol.put(i, ts);
     591      }
     592    }
     593  }
     594  return out;
     595}
     596
     597CountedPtr< Scantable > STMath::arrayOperateRow( const CountedPtr< Scantable >& in,
     598                                                 const std::vector<float> val,
     599                                                 const std::string& mode,
     600                                                 bool tsys )
     601{
     602  // conformity of SPECTRA and TSYS
     603  if ( tsys ) {
     604    TableIterator titer(in->table(), "IFNO");
     605    while ( !titer.pastEnd() ) {
     606      ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
     607      ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
     608      Array<Float> spec = specCol.getColumn() ;
     609      Array<Float> ts = tsysCol.getColumn() ;
     610      if ( !spec.conform( ts ) ) {
     611        throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
     612      }
     613      titer.next() ;
     614    }
     615  }
     616
     617  // check if vector size is equal to nrow
     618  Vector<Float> fact( val ) ;
     619  if ( fact.nelements() != in->nrow() ) {
     620    throw( AipsError("Vector size must be same as number of row.") ) ;
     621  }
     622
     623  // check divided by zero
     624  if ( ( mode == "DIV" ) && anyEQ( fact, (float)0.0 ) ) {
     625    throw( AipsError("Divided by zero is not recommended." ) ) ;
     626  }
     627
     628  CountedPtr< Scantable > out = getScantable(in, false);
     629  Table& tab = out->table();
     630  ArrayColumn<Float> specCol(tab,"SPECTRA");
     631  ArrayColumn<Float> tsysCol(tab,"TSYS");
     632  if (mode == "DIV") fact = (float)1.0 / fact;
     633  if (mode == "SUB") fact *= (float)-1.0 ;
     634  for (uInt i=0; i<tab.nrow(); ++i) {
     635    Vector<Float> spec;
     636    Vector<Float> ts;
     637    specCol.get(i, spec);
     638    tsysCol.get(i, ts);
     639    if (mode == "MUL" || mode == "DIV") {
     640      spec *= fact[i];
     641      specCol.put(i, spec);
     642      if ( tsys ) {
     643        ts *= fact[i];
     644        tsysCol.put(i, ts);
     645      }
     646    } else if ( mode == "ADD"  || mode == "SUB") {
     647      spec += fact[i];
     648      specCol.put(i, spec);
     649      if ( tsys ) {
     650        ts += fact[i];
     651        tsysCol.put(i, ts);
     652      }
     653    }
     654  }
     655  return out;
     656}
     657
     658CountedPtr< Scantable > STMath::array2dOperate( const CountedPtr< Scantable >& in,
     659                                                const std::vector< std::vector<float> > val,
     660                                                const std::string& mode,
     661                                                bool tsys )
     662{
     663  // conformity of SPECTRA and TSYS
     664  if ( tsys ) {
     665    TableIterator titer(in->table(), "IFNO");
     666    while ( !titer.pastEnd() ) {
     667      ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
     668      ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
     669      Array<Float> spec = specCol.getColumn() ;
     670      Array<Float> ts = tsysCol.getColumn() ;
     671      if ( !spec.conform( ts ) ) {
     672        throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
     673      }
     674      titer.next() ;
     675    }
     676  }
     677
     678  // some checks
     679  vector<uInt> nchans;
     680  for ( uInt i = 0 ; i < in->nrow() ; i++ ) {
     681    nchans.push_back( (in->getSpectrum( i )).size() ) ;
     682  }
     683  //Vector<uInt> mchans( nchans ) ;
     684  vector< Vector<Float> > facts ;
     685  for ( uInt i = 0 ; i < nchans.size() ; i++ ) {
     686    Vector<Float> tmp( val[i] ) ;
     687    // check divided by zero
     688    if ( ( mode == "DIV" ) && anyEQ( tmp, (float)0.0 ) ) {
     689      throw( AipsError("Divided by zero is not recommended." ) ) ;
     690    }
     691    // conformity check
     692    if ( tmp.nelements() != nchans[i] ) {
     693      stringstream ss ;
     694      ss << "Row " << i << ": Vector size must be same as number of channel." ;
     695      throw( AipsError( ss.str() ) ) ;
     696    }
     697    facts.push_back( tmp ) ;
     698  }
     699
     700
     701  CountedPtr< Scantable > out = getScantable(in, false);
     702  Table& tab = out->table();
     703  ArrayColumn<Float> specCol(tab,"SPECTRA");
     704  ArrayColumn<Float> tsysCol(tab,"TSYS");
     705  for (uInt i=0; i<tab.nrow(); ++i) {
     706    Vector<Float> fact = facts[i] ;
     707    Vector<Float> spec;
     708    Vector<Float> ts;
     709    specCol.get(i, spec);
     710    tsysCol.get(i, ts);
     711    if (mode == "MUL" || mode == "DIV") {
     712      if (mode == "DIV") fact = (float)1.0 / fact;
     713      spec *= fact;
     714      specCol.put(i, spec);
     715      if ( tsys ) {
     716        ts *= fact;
     717        tsysCol.put(i, ts);
     718      }
     719    } else if ( mode == "ADD"  || mode == "SUB") {
     720      if (mode == "SUB") fact *= (float)-1.0 ;
     721      spec += fact;
     722      specCol.put(i, spec);
     723      if ( tsys ) {
     724        ts += fact;
    499725        tsysCol.put(i, ts);
    500726      }

  • branches/alma/src/STMath.h

    r1673 r1677  
    117117  casa::CountedPtr<Scantable>
    118118    unaryOperate( const casa::CountedPtr<Scantable>& in, float val,
     119                  const std::string& mode, bool tsys=false );
     120
     121  // array operation
     122  casa::CountedPtr<Scantable>
     123    arrayOperate( const casa::CountedPtr<Scantable>& in,
     124                  const std::vector<float> val,
     125                  const std::string& mode, bool tsys=false,
     126                  const std::string& op="channel" );
     127
     128  // channel operation
     129  casa::CountedPtr<Scantable>
     130    arrayOperateChannel( const casa::CountedPtr<Scantable>& in,
     131                         const std::vector<float> val,
     132                         const std::string& mode, bool tsys=false );
     133
     134  // row operation
     135  casa::CountedPtr<Scantable>
     136    arrayOperateRow( const casa::CountedPtr<Scantable>& in,
     137                     const std::vector<float> val,
     138                     const std::string& mode, bool tsys=false );
     139
     140  // 2d array operation
     141  casa::CountedPtr<Scantable>
     142    array2dOperate( const casa::CountedPtr<Scantable>& in,
     143                  const std::vector< std::vector<float> > val,
    119144                  const std::string& mode, bool tsys=false );
    120145

  • branches/alma/src/STMathWrapper.h

    r1673 r1677  
    7373  { return ScantableWrapper(STMath::unaryOperate(in.getCP(), val, mode, tsys)); }
    7474
     75  ScantableWrapper arrayOperate( const ScantableWrapper& in,
     76                                 const std::vector<float> val,
     77                                 const std::string& mode, bool tsys=false,
     78                                 const std::string& op="channel" )
     79  { return ScantableWrapper(STMath::arrayOperate(in.getCP(), val, mode, tsys, op)); }
     80
     81  ScantableWrapper array2dOperate( const ScantableWrapper& in,
     82                                   const std::vector< std::vector<float> > val,
     83                                   const std::string& mode, bool tsys=false )
     84  { return ScantableWrapper(STMath::array2dOperate(in.getCP(), val, mode, tsys)); }
     85
    7586  ScantableWrapper binaryOperate( const ScantableWrapper& left,
    7687                                  const ScantableWrapper& right,

  • branches/alma/src/python_STMath.cpp

    r1673 r1677  
    4949        .def("_averagebeams", &STMathWrapper::averageBeams)
    5050        .def("_unaryop", &STMathWrapper::unaryOperate)
     51        .def("_arrayop", &STMathWrapper::arrayOperate)
     52        //.def("_array2dop", &STMathWrapper::array2dOperate)
    5153        .def("_binaryop", &STMathWrapper::binaryOperate)
    5254        .def("_auto_quotient", &STMathWrapper::autoQuotient)
Note: See TracChangeset for help on using the changeset viewer.