Changeset 165 for trunk


Ignore:
Timestamp:
12/27/04 20:31:36 (20 years ago)
Author:
kil064
Message:

Reimplement 'average_pol' with both insitu and 'outsit' versions

Location:
trunk/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/SDMath.cc

    r163 r165  
    448448}
    449449
    450 
    451 CountedPtr<SDMemTable>
    452 SDMath::averagePol(const CountedPtr<SDMemTable>& in,
    453                    const Vector<Bool>& mask)
    454 //
    455 // Average all polarizations together, weighted by variance
    456 //
    457 {
    458 //   WeightType wtType = NONE;
    459 //   convertWeightString (wtType, weight);
    460 
    461    const uInt nRows = in->nRow();
    462    const uInt polAxis = 2;                     // Polarization axis
    463    const uInt chanAxis = 3;                    // Spectrum axis
    464 
    465 // Create output Table and reshape number of polarizations
    466 
    467   Bool clear=True;
    468   SDMemTable* pTabOut = new SDMemTable(*in, clear);
    469   SDHeader header = pTabOut->getSDHeader();
    470   header.npol = 1;
    471   pTabOut->putSDHeader(header);
    472 
    473 // Shape of input and output data
    474 
    475   const IPosition shapeIn = in->rowAsMaskedArray(0).shape();
    476   IPosition shapeOut(shapeIn);
    477   shapeOut(polAxis) = 1;                          // Average all polarizations
    478 //
    479   const uInt nChan = shapeIn(chanAxis);
    480   const IPosition vecShapeOut(4,1,1,1,nChan);     // A multi-dim form of a Vector shape
    481   IPosition start(4), end(4);
    482 
    483 // Output arrays
    484 
    485   Array<Float> outData(shapeOut, 0.0);
    486   Array<Bool> outMask(shapeOut, True);
    487   const IPosition axes(2, 2, 3);              // pol-channel plane
    488 //
    489   const Bool useMask = (mask.nelements() == shapeIn(chanAxis));
    490 
    491 // Loop over rows
    492 
    493    for (uInt iRow=0; iRow<nRows; iRow++) {
    494 
    495 // Get data for this row
    496 
    497       MaskedArray<Float> marr(in->rowAsMaskedArray(iRow));
    498       Array<Float>& arr = marr.getRWArray();
    499       const Array<Bool>& barr = marr.getMask();
    500 
    501 // Make iterators to iterate by pol-channel planes
    502 
    503       ReadOnlyArrayIterator<Float> itDataPlane(arr, axes);
    504       ReadOnlyArrayIterator<Bool> itMaskPlane(barr, axes);
    505 
    506 // Accumulations
    507 
    508       Float fac = 1.0;
    509       Vector<Float> vecSum(nChan,0.0);
    510 
    511 // Iterate through data by pol-channel planes
    512 
    513       while (!itDataPlane.pastEnd()) {
    514 
    515 // Iterate through plane by polarization  and accumulate Vectors
    516 
    517         Vector<Float> t1(nChan); t1 = 0.0;
    518         Vector<Bool> t2(nChan); t2 = True;
    519         MaskedArray<Float> vecSum(t1,t2);
    520         Float varSum = 0.0;
    521         {
    522            ReadOnlyVectorIterator<Float> itDataVec(itDataPlane.array(), 1);
    523            ReadOnlyVectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
    524            while (!itDataVec.pastEnd()) {     
    525 
    526 // Create MA of data & mask (optionally including OTF mask) and  get variance
    527 
    528               if (useMask) {
    529                  const MaskedArray<Float> spec(itDataVec.vector(),mask&&itMaskVec.vector());
    530                  fac = 1.0 / variance(spec);
    531               } else {
    532                  const MaskedArray<Float> spec(itDataVec.vector(),itMaskVec.vector());
    533                  fac = 1.0 / variance(spec);
    534               }
    535 
    536 // Normalize spectrum (without OTF mask) and accumulate
    537 
    538               const MaskedArray<Float> spec(fac*itDataVec.vector(), itMaskVec.vector());
    539               vecSum += spec;
    540               varSum += fac;
    541 
    542 // Next
    543 
    544               itDataVec.next();
    545               itMaskVec.next();
    546            }
    547         }
    548 
    549 // Normalize summed spectrum
    550 
    551         vecSum /= varSum;
    552 
    553 // FInd position in input data array.  We are iterating by pol-channel
    554 // plane so all that will change is beam and IF and that's what we want.
    555 
    556         IPosition pos = itDataPlane.pos();
    557 
    558 // Write out data. This is a bit messy. We have to reform the Vector
    559 // accumulator into an Array of shape (1,1,1,nChan)
    560 
    561         start = pos;
    562         end = pos;
    563         end(chanAxis) = nChan-1;
    564         outData(start,end) = vecSum.getArray().reform(vecShapeOut);
    565         outMask(start,end) = vecSum.getMask().reform(vecShapeOut);
    566 
    567 // Step to next beam/IF combination
    568 
    569         itDataPlane.next();
    570         itMaskPlane.next();
    571       }
    572 
    573 // Generate output container and write it to output table
    574 
    575       SDContainer sc = in->getSDContainer();
    576       sc.resize(shapeOut);
    577 //
    578       putDataInSDC (sc, outData, outMask);
    579       pTabOut->putSDContainer(sc);
    580    }
    581 //
    582   return CountedPtr<SDMemTable>(pTabOut);
    583 }
    584 
     450void SDMath::averagePolInSitu (SDMemTable* pIn, const Vector<Bool>& mask)
     451{
     452  SDMemTable* pOut = localAveragePol (*pIn, mask);
     453  *pIn = *pOut;
     454   delete pOut;
     455} 
     456
     457
     458CountedPtr<SDMemTable> SDMath::averagePol (const CountedPtr<SDMemTable>& in,
     459                                           const Vector<Bool>& mask)
     460{
     461  return CountedPtr<SDMemTable>(localAveragePol(*in, mask));
     462}
    585463
    586464CountedPtr<SDMemTable> SDMath::bin(const CountedPtr<SDMemTable>& in,
    … …  
    949827    sc.putFlags(outflags);
    950828}
     829
     830
     831SDMemTable* SDMath::localAveragePol(const SDMemTable& in, const Vector<Bool>& mask)
     832//
     833// Average all polarizations together, weighted by variance
     834//
     835{
     836//   WeightType wtType = NONE;
     837//   convertWeightString (wtType, weight);
     838
     839   const uInt nRows = in.nRow();
     840   const uInt polAxis = 2;                     // Polarization axis
     841   const uInt chanAxis = 3;                    // Spectrum axis
     842
     843// Create output Table and reshape number of polarizations
     844
     845  Bool clear=True;
     846  SDMemTable* pTabOut = new SDMemTable(in, clear);
     847  SDHeader header = pTabOut->getSDHeader();
     848  header.npol = 1;
     849  pTabOut->putSDHeader(header);
     850
     851// Shape of input and output data
     852
     853  const IPosition& shapeIn = in.rowAsMaskedArray(0u, False).shape();
     854  IPosition shapeOut(shapeIn);
     855  shapeOut(polAxis) = 1;                          // Average all polarizations
     856//
     857  const uInt nChan = shapeIn(chanAxis);
     858  const IPosition vecShapeOut(4,1,1,1,nChan);     // A multi-dim form of a Vector shape
     859  IPosition start(4), end(4);
     860
     861// Output arrays
     862
     863  Array<Float> outData(shapeOut, 0.0);
     864  Array<Bool> outMask(shapeOut, True);
     865  const IPosition axes(2, 2, 3);              // pol-channel plane
     866//
     867  const Bool useMask = (mask.nelements() == shapeIn(chanAxis));
     868cerr << "nEl=" << mask.nelements() << endl;
     869cerr << "useMask=" << useMask << endl;
     870
     871// Loop over rows
     872
     873   for (uInt iRow=0; iRow<nRows; iRow++) {
     874
     875// Get data for this row
     876
     877      MaskedArray<Float> marr(in.rowAsMaskedArray(iRow));
     878      Array<Float>& arr = marr.getRWArray();
     879      const Array<Bool>& barr = marr.getMask();
     880
     881// Make iterators to iterate by pol-channel planes
     882
     883      ReadOnlyArrayIterator<Float> itDataPlane(arr, axes);
     884      ReadOnlyArrayIterator<Bool> itMaskPlane(barr, axes);
     885
     886// Accumulations
     887
     888      Float fac = 1.0;
     889      Vector<Float> vecSum(nChan,0.0);
     890
     891// Iterate through data by pol-channel planes
     892
     893      while (!itDataPlane.pastEnd()) {
     894
     895// Iterate through plane by polarization  and accumulate Vectors
     896
     897        Vector<Float> t1(nChan); t1 = 0.0;
     898        Vector<Bool> t2(nChan); t2 = True;
     899        MaskedArray<Float> vecSum(t1,t2);
     900        Float varSum = 0.0;
     901        {
     902           ReadOnlyVectorIterator<Float> itDataVec(itDataPlane.array(), 1);
     903           ReadOnlyVectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
     904           while (!itDataVec.pastEnd()) {     
     905
     906// Create MA of data & mask (optionally including OTF mask) and  get variance
     907
     908              if (useMask) {
     909                 const MaskedArray<Float> spec(itDataVec.vector(),mask&&itMaskVec.vector());
     910                 fac = 1.0 / variance(spec);
     911              } else {
     912                 const MaskedArray<Float> spec(itDataVec.vector(),itMaskVec.vector());
     913                 fac = 1.0 / variance(spec);
     914              }
     915
     916// Normalize spectrum (without OTF mask) and accumulate
     917
     918              const MaskedArray<Float> spec(fac*itDataVec.vector(), itMaskVec.vector());
     919              vecSum += spec;
     920              varSum += fac;
     921
     922// Next
     923
     924              itDataVec.next();
     925              itMaskVec.next();
     926           }
     927        }
     928
     929// Normalize summed spectrum
     930
     931        vecSum /= varSum;
     932
     933// FInd position in input data array.  We are iterating by pol-channel
     934// plane so all that will change is beam and IF and that's what we want.
     935
     936        IPosition pos = itDataPlane.pos();
     937
     938// Write out data. This is a bit messy. We have to reform the Vector
     939// accumulator into an Array of shape (1,1,1,nChan)
     940
     941        start = pos;
     942        end = pos;
     943        end(chanAxis) = nChan-1;
     944        outData(start,end) = vecSum.getArray().reform(vecShapeOut);
     945        outMask(start,end) = vecSum.getMask().reform(vecShapeOut);
     946
     947// Step to next beam/IF combination
     948
     949        itDataPlane.next();
     950        itMaskPlane.next();
     951      }
     952
     953// Generate output container and write it to output table
     954
     955      SDContainer sc = in.getSDContainer();
     956      sc.resize(shapeOut);
     957//
     958      putDataInSDC (sc, outData, outMask);
     959      pTabOut->putSDContainer(sc);
     960   }
     961//
     962  return pTabOut;
     963}

  • trunk/src/SDMath.h

    r162 r165  
    6767
    6868  casa::CountedPtr<SDMemTable> bin(const casa::CountedPtr<SDMemTable>& in,
    69                              casa::Int width);
     69                                   casa::Int width);
    7070
    7171// Average in time
    7272
    73   casa::CountedPtr<SDMemTable>
    74   average (const casa::Block<casa::CountedPtr<SDMemTable> >& in,
    75            const casa::Vector<casa::Bool>& mask,
    76            bool scanAverage, const std::string& weightStr);
     73  casa::CountedPtr<SDMemTable>  average (const casa::Block<casa::CountedPtr<SDMemTable> >& in,
     74                                         const casa::Vector<casa::Bool>& mask,
     75                                         bool scanAverage, const std::string& weightStr);
    7776
    7877// Average polarizations
    7978
    80   casa::CountedPtr<SDMemTable>
    81   averagePol(const casa::CountedPtr<SDMemTable>& in, const casa::Vector<casa::Bool>& mask);
     79  void averagePolInSitu (SDMemTable* in, const casa::Vector<casa::Bool>& mask);
     80  casa::CountedPtr<SDMemTable> averagePol(const casa::CountedPtr<SDMemTable>& in,
     81                                          const casa::Vector<casa::Bool>& mask);
    8282
    8383// Statistics
    … …  
    135135
    136136   void convertWeightString (WeightType& wt, const std::string& weightStr);
     137
     138// Function for simple mathematical operations.  what=0 (mul) or 1 (add)
     139
     140  SDMemTable* localOperate (const SDMemTable& in, casa::Float offset,
     141                            casa::Bool doAll, casa::uInt what);
     142
     143// Function to average polarizations
     144
     145  SDMemTable* localAveragePol(const SDMemTable& in, const casa::Vector<casa::Bool>& mask);
    137146};
    138147
    … …  
    140149
    141150#endif
    142 
    143 
    144 
    145 
    146 

  • trunk/src/SDMathWrapper.h

    r151 r165  
    9898// Average polarizations
    9999
    100   SDMemTableWrapper averagePol(const SDMemTableWrapper& in,
    101                                const std::vector<bool>& mask) {
    102     return SDMath::averagePol(in.getCP(), mask);
     100  void averagePolInSitu (SDMemTableWrapper& in,  const std::vector<bool>& mask)
     101  {
     102    SDMemTable* sdmt = in.getPtr();
     103    SDMath::averagePolInSitu(in.getPtr(), mask);
     104  }
     105  SDMemTableWrapper averagePol (const SDMemTableWrapper& in, const std::vector<bool>& mask)
     106  {
     107    return SDMemTableWrapper(SDMath::averagePol(in.getCP(), mask));
    103108  }
    104109

  • trunk/src/python_SDMath.cc

    r153 r165  
    6969      def("average", &SDMathWrapper::average);
    7070      def("averagepol", &SDMathWrapper::averagePol);
     71      def("averagepol_insitu", &SDMathWrapper::averagePolInSitu);
    7172      def("bin", &SDMathWrapper::bin);
    7273      def("stats", &SDMathWrapper::statistic);
Note: See TracChangeset for help on using the changeset viewer.