Changeset 146 for trunk


Ignore:
Timestamp:
12/26/04 14:56:06 (20 years ago)
Author:
kil064
Message:

rework 'multiply' and 'multiplyInSitu' to use one common
'localMultiply' function behind the scenes.

Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/SDMath.cc

    r144 r146  
    268268        }
    269269
    270 // Accumulation step. First get data and deconstruct
    271 
    272         MaskedArray<Float> dataIn(in[iTab]->rowAsMaskedArray(iRow));
    273         Array<Float>& valuesIn = dataIn.getRWArray();           // writable reference
    274         const Array<Bool>& maskIn = dataIn.getMask();          // RO reference
    275 //
    276         if (wtType==NONE) {
    277            const MaskedArray<Float> n(nInc,dataIn.getMask());
    278            nPts += n;                               // Only accumulates where mask==T
    279         } else if (wtType==VAR) {
    280 
    281 // We are going to average the data, weighted by the noise for each pol, beam and IF.
    282 // So therefore we need to iterate through by spectrum (axis 3)
    283 
    284            VectorIterator<Float> itData(valuesIn, axis);
    285            ReadOnlyVectorIterator<Bool> itMask(maskIn, axis);
    286            while (!itData.pastEnd()) {
    287 
    288 // Make MaskedArray of Vector, optionally apply OTF mask, and find scaling factor
    289 
    290              if (useMask) {
    291                 MaskedArray<Float> tmp(itData.vector(),mask&&itMask.vector());
    292                 fac = 1.0/variance(tmp);
    293              } else {
    294                 MaskedArray<Float> tmp(itData.vector(),itMask.vector());
    295                 fac = 1.0/variance(tmp);
    296              }
    297 
    298 // Scale data
    299 
    300              itData.vector() *= fac;     // Writes back into 'dataIn'
    301 //
    302 // Accumulate variance per if/pol/beam averaged over spectrum
    303 // This method to get pos2 from itData.pos() is only valid
    304 // because the spectral axis is the last one (so we can just
    305 // copy the first nAXesSub positions out)
    306 
    307              pos2 = itData.pos().getFirst(nAxesSub);
    308              sumSq(pos2) += fac;
    309 //
    310              itData.next();
    311              itMask.next();
    312            }
    313         } else if (wtType==TSYS) {
    314         }
    315 
    316 // Accumulate sum of (possibly scaled) data
    317 
    318        sum += dataIn;
    319 
    320 // Accumulate Tsys, time, and interval
    321 
    322        tSysSum += tSys;
    323        timeSum += time;
    324        intSum += interval;
    325 
    326 // Number of rows in accumulation
    327 
    328        nAccum += 1;
     270// Accumulate
     271
     272        accumulate (timeSum, intSum, nAccum, sum, sumSq, nPts, tSysSum,
     273                    tSys, nInc, mask, time, interval, in, iTab, iRow, axis,
     274                    nAxesSub, useMask, wtType);
     275//
    329276       oldSourceName = sourceName;
    330277       oldFreqID = freqID;
    … …  
    433380}
    434381
    435 void SDMath::multiplyInSitu(SDMemTable* in, Float factor) {
    436   SDMemTable* sdmt = new SDMemTable(*in);
    437   Table t = sdmt->table();
    438   ArrayColumn<Float> spec(t,"SPECTRA"); 
    439   for (uInt i=0; i < t.nrow(); i++) {
    440     MaskedArray<Float> marr(sdmt->rowAsMaskedArray(i));
    441     marr *= factor;
    442     spec.put(i, marr.getArray());
    443   }
    444   in = sdmt;
    445   delete sdmt;sdmt=0;
    446 }
     382
     383
     384void SDMath::multiplyInSitu(SDMemTable* pIn, Float factor)
     385{
     386  SDMemTable* pOut = localMultiply (*pIn, factor);
     387  *pIn = *pOut;
     388   delete pOut;
     389} 
     390
    447391
    448392CountedPtr<SDMemTable>
    449393SDMath::multiply(const CountedPtr<SDMemTable>& in, Float factor)
    450 //
    451 // Multiply values by factor
    452 //
    453 {
    454   SDMemTable* sdmt = new SDMemTable(*in);
    455   Table t = sdmt->table();
    456   ArrayColumn<Float> spec(t,"SPECTRA");
    457 
    458   for (uInt i=0; i < t.nrow(); i++) {
    459     MaskedArray<Float> marr(sdmt->rowAsMaskedArray(i));
    460     marr *= factor;
    461     spec.put(i, marr.getArray());
    462   }
    463   return CountedPtr<SDMemTable>(sdmt);
     394{
     395  return CountedPtr<SDMemTable>(localMultiply(*in,factor));
    464396}
    465397
    … …  
    470402//
    471403{
    472   SDMemTable* sdmt = new SDMemTable(*in);
    473 
     404  SDMemTable* sdmt = new SDMemTable(*in,False);
    474405  Table t = sdmt->table();
    475406  ArrayColumn<Float> spec(t,"SPECTRA");
    … …  
    765696}
    766697
     698
     699
     700// 'private' functions
     701
    767702void SDMath::fillSDC (SDContainer& sc,
    768703                      const Array<Bool>& mask,
    … …  
    823758}
    824759
     760
     761void SDMath::accumulate (Double& timeSum, Double& intSum, Int& nAccum,
     762                         MaskedArray<Float>& sum, Array<Float>& sumSq,
     763                         Array<Float>& nPts, Array<Float>& tSysSum,
     764                         const Array<Float>& tSys, const Array<Float>& nInc,
     765                         const Vector<Bool>& mask, Double time, Double interval,
     766                         const Block<CountedPtr<SDMemTable> >& in,
     767                         uInt iTab, uInt iRow, uInt axis,
     768                         uInt nAxesSub, Bool useMask,
     769                         weightType wtType)
     770{
     771
     772// Get data
     773
     774   MaskedArray<Float> dataIn(in[iTab]->rowAsMaskedArray(iRow));
     775   Array<Float>& valuesIn = dataIn.getRWArray();           // writable reference
     776   const Array<Bool>& maskIn = dataIn.getMask();          // RO reference
     777//
     778   if (wtType==NONE) {
     779      const MaskedArray<Float> n(nInc,dataIn.getMask());
     780      nPts += n;                               // Only accumulates where mask==T
     781   } else if (wtType==VAR) {
     782
     783// We are going to average the data, weighted by the noise for each pol, beam and IF.
     784// So therefore we need to iterate through by spectrum (axis 3)
     785
     786      VectorIterator<Float> itData(valuesIn, axis);
     787      ReadOnlyVectorIterator<Bool> itMask(maskIn, axis);
     788      Float fac = 1.0;
     789      IPosition pos(nAxesSub,0); 
     790//
     791      while (!itData.pastEnd()) {
     792
     793// Make MaskedArray of Vector, optionally apply OTF mask, and find scaling factor
     794
     795        if (useMask) {
     796           MaskedArray<Float> tmp(itData.vector(),mask&&itMask.vector());
     797           fac = 1.0/variance(tmp);
     798        } else {
     799           MaskedArray<Float> tmp(itData.vector(),itMask.vector());
     800           fac = 1.0/variance(tmp);
     801        }
     802
     803// Scale data
     804
     805        itData.vector() *= fac;     // Writes back into 'dataIn'
     806//
     807// Accumulate variance per if/pol/beam averaged over spectrum
     808// This method to get pos2 from itData.pos() is only valid
     809// because the spectral axis is the last one (so we can just
     810// copy the first nAXesSub positions out)
     811
     812        pos = itData.pos().getFirst(nAxesSub);
     813        sumSq(pos) += fac;
     814//
     815        itData.next();
     816        itMask.next();
     817      }
     818   } else if (wtType==TSYS) {
     819   }
     820
     821// Accumulate sum of (possibly scaled) data
     822
     823   sum += dataIn;
     824
     825// Accumulate Tsys, time, and interval
     826
     827   tSysSum += tSys;
     828   timeSum += time;
     829   intSum += interval;
     830   nAccum += 1;
     831}
     832
     833SDMemTable* SDMath::localMultiply (const SDMemTable& in, Float factor)
     834{
     835  SDMemTable* pOut = new SDMemTable(in,False);
     836  const Table& tOut = pOut->table();
     837  ArrayColumn<Float> spec(tOut,"SPECTRA"); 
     838//
     839  for (uInt i=0; i < tOut.nrow(); i++) {
     840    MaskedArray<Float> marr(pOut->rowAsMaskedArray(i));
     841    marr *= factor;
     842    spec.put(i, marr.getArray());
     843  }
     844   return pOut;
     845}
     846
     847

  • trunk/src/SDMath.h

    r144 r146  
    4545  casa::CountedPtr<SDMemTable> quotient(const casa::CountedPtr<SDMemTable>& on,
    4646                                         const casa::CountedPtr<SDMemTable>& off);
    47 
     47//
    4848  void multiplyInSitu(SDMemTable* in, casa::Float factor);
    49 
    5049  casa::CountedPtr<SDMemTable> multiply(const casa::CountedPtr<SDMemTable>& in,
    51                                   casa::Float factor);
    52 
     50                                        casa::Float factor);
     51//
    5352  casa::CountedPtr<SDMemTable> add(const casa::CountedPtr<SDMemTable>& in,
    5453                             casa::Float offset);
    … …  
    7473  enum weightType {NONE,VAR,TSYS};
    7574
     75  void accumulate (casa::Double& timeSum, casa::Double& intSum, casa::Int& nAccum,
     76                   casa::MaskedArray<casa::Float>& sum, casa::Array<casa::Float>& sumSq,
     77                   casa::Array<casa::Float>& nPts, casa::Array<casa::Float>& tSysSum,
     78                   const casa::Array<casa::Float>& tSys,  const casa::Array<casa::Float>& nInc,
     79                   const casa::Vector<casa::Bool>& mask, casa::Double time, casa::Double interval,
     80                   const casa::Block<casa::CountedPtr<SDMemTable> >& in,
     81                   casa::uInt iTab, casa::uInt iRow, casa::uInt axis, casa::uInt nAxesSub,
     82                   casa::Bool useMask, weightType wtType);
     83
    7684  void fillSDC (SDContainer& sc, const casa::Array<casa::Bool>& mask,
    7785                const casa::Array<casa::Float>& data,
    … …  
    8088                casa::Double interval, const casa::String& sourceName,
    8189                const casa::Vector<casa::uInt>& freqID);
    82    void normalize (casa::MaskedArray<casa::Float>& data,
    83                    const casa::Array<casa::Float>& sumSq,
    84                    const casa::Array<casa::Float>& nPts,
    85                    weightType wtType, casa::Int axis, casa::Int nAxes);
    8690
     91  SDMemTable* localMultiply (const SDMemTable& in, casa::Float factor);
    8792
     93  void normalize (casa::MaskedArray<casa::Float>& data,
     94                  const casa::Array<casa::Float>& sumSq,
     95                  const casa::Array<casa::Float>& nPts,
     96                  weightType wtType, casa::Int axis, casa::Int nAxes);
    8897};
    8998
    … …  
    91100
    92101#endif
     102
     103
     104
     105
     106
Note: See TracChangeset for help on using the changeset viewer.