Changeset 518


Ignore:
Timestamp:
02/28/05 21:48:50 (20 years ago)
Author:
kil064
Message:

add Tsys and Tint weightring

Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/SDMath.cc

    r502 r518  
    150150 
    151151  WeightType wtType = NONE;
    152   convertWeightString(wtType, weightStr);
     152  convertWeightString(wtType, weightStr, True);
    153153
    154154// Create output Table by cloning from the first table
     
    219219// To get the right shape for the Tsys accumulator we need to
    220220// access a column from the first table.  The shape of this
    221 // array must not change
    222 
    223   Array<Float> tSysSum;
     221// array must not change.  Note however that since the TSysSqSum
     222// array is used in a normalization process, and that I ignore the
     223// channel axis replication of values for now, it loses a dimension
     224
     225  Array<Float> tSysSum, tSysSqSum;
    224226  {
    225227    const Table& tabIn = in[0]->table();
    226228    tSysCol.attach(tabIn,"TSYS");
    227229    tSysSum.resize(tSysCol.shape(0));
     230//
     231    tSysSqSum.resize(shp2);
    228232  }
    229233  tSysSum =0.0;
     234  tSysSqSum = 0.0;
    230235  Array<Float> tSys;
    231236
     
    306311// Normalize data in 'sum' accumulation array according to weighting scheme
    307312
    308            normalize(sum, sumSq, nPts, wtType, asap::ChanAxis, nAxesSub);
     313           normalize(sum, sumSq, tSysSqSum, nPts, intSum, wtType, asap::ChanAxis, nAxesSub);
    309314
    310315// Get ScanContainer for the first row of this averaged Scan
     
    331336//
    332337           tSysSum =0.0;
     338           tSysSqSum =0.0;
    333339           timeSum = 0.0;
    334340           intSum = 0.0;
     
    348354// Accumulate
    349355
    350         accumulate(timeSum, intSum, nAccum, sum, sumSq, nPts, tSysSum,
    351                     tSys, nInc, mask, time, interval, in, iTab, iRow, asap::ChanAxis,
    352                     nAxesSub, useMask, wtType);
     356        accumulate(timeSum, intSum, nAccum, sum, sumSq, nPts, tSysSum, tSysSqSum,
     357                   tSys, nInc, mask, time, interval, in, iTab, iRow, asap::ChanAxis,
     358                   nAxesSub, useMask, wtType);
    353359//
    354360       oldSourceName = sourceName;
     
    364370// Normalize data in 'sum' accumulation array according to weighting scheme
    365371
    366   normalize(sum, sumSq, nPts, wtType, asap::ChanAxis, nAxesSub);
     372  normalize(sum, sumSq, tSysSqSum, nPts, intSum, wtType, asap::ChanAxis, nAxesSub);
    367373
    368374// Create and fill container.  The container we clone will be from
     
    807813{
    808814   WeightType wtType = NONE;
    809    convertWeightString(wtType, weightStr);
     815   convertWeightString(wtType, weightStr, False);
    810816
    811817   const uInt nRows = in.nRow();
     
    12321238//
    12331239// phase in degrees
    1234 // Applies to all Beams and IFs
    1235 // Might want to optionally select on Beam/IF
     1240// assumes linear correlations
    12361241//
    12371242{
    12381243   if (in.nPol() != 4) {
    12391244      throw(AipsError("You must have 4 polarizations to run this function"));
     1245   }
     1246//
     1247   SDHeader sh = in.getSDHeader();
     1248   Instrument inst = SDAttr::convertInstrument (sh.antennaname, False);
     1249   SDAttr sdAtt;
     1250   if (sdAtt.feedPolType(inst) != LINEAR) {
     1251      throw(AipsError("Only linear polarizations are supported"));
    12401252   }
    12411253//   
     
    12841296//
    12851297// phase in degrees
    1286 // Applies to all Beams and IFs
    1287 // Might want to optionally select on Beam/IF
     1298// assumes linear correlations
    12881299//
    12891300{
    12901301   if (in.nPol() != 4) {
    12911302      throw(AipsError("You must have 4 polarizations to run this function"));
     1303   }
     1304//
     1305   SDHeader sh = in.getSDHeader();
     1306   Instrument inst = SDAttr::convertInstrument (sh.antennaname, False);
     1307   SDAttr sdAtt;
     1308   if (sdAtt.feedPolType(inst) != LINEAR) {
     1309      throw(AipsError("Only linear polarizations are supported"));
    12921310   }
    12931311//   
     
    17001718}
    17011719
     1720void SDMath::accumulate(Double& timeSum, Double& intSum, Int& nAccum,
     1721                        MaskedArray<Float>& sum, Array<Float>& sumSq,
     1722                        Array<Float>& nPts, Array<Float>& tSysSum,
     1723                        Array<Float>& tSysSqSum,
     1724                        const Array<Float>& tSys, const Array<Float>& nInc,
     1725                        const Vector<Bool>& mask, Double time, Double interval,
     1726                        const Block<CountedPtr<SDMemTable> >& in,
     1727                        uInt iTab, uInt iRow, uInt axis,
     1728                        uInt nAxesSub, Bool useMask,
     1729                        WeightType wtType) const
     1730{
     1731
     1732// Get data
     1733
     1734   MaskedArray<Float> dataIn(in[iTab]->rowAsMaskedArray(iRow));
     1735   Array<Float>& valuesIn = dataIn.getRWArray();           // writable reference
     1736   const Array<Bool>& maskIn = dataIn.getMask();          // RO reference
     1737//
     1738   if (wtType==NONE) {
     1739      const MaskedArray<Float> n(nInc,dataIn.getMask());
     1740      nPts += n;                               // Only accumulates where mask==T
     1741   } else if (wtType==TINT) {
     1742
     1743// We are weighting the data by integration time.
     1744
     1745     valuesIn *= Float(interval);
     1746
     1747   } else if (wtType==VAR) {
     1748
     1749// We are going to average the data, weighted by the noise for each pol, beam and IF.
     1750// So therefore we need to iterate through by spectrum (axis 3)
     1751
     1752      VectorIterator<Float> itData(valuesIn, axis);
     1753      ReadOnlyVectorIterator<Bool> itMask(maskIn, axis);
     1754      Float fac = 1.0;
     1755      IPosition pos(nAxesSub,0); 
     1756//
     1757      while (!itData.pastEnd()) {
     1758
     1759// Make MaskedArray of Vector, optionally apply OTF mask, and find scaling factor
     1760
     1761         if (useMask) {
     1762            MaskedArray<Float> tmp(itData.vector(),mask&&itMask.vector());
     1763            fac = 1.0/variance(tmp);
     1764         } else {
     1765            MaskedArray<Float> tmp(itData.vector(),itMask.vector());
     1766            fac = 1.0/variance(tmp);
     1767         }
     1768
     1769// Scale data
     1770
     1771         itData.vector() *= fac;     // Writes back into 'dataIn'
     1772//
     1773// Accumulate variance per if/pol/beam averaged over spectrum
     1774// This method to get pos2 from itData.pos() is only valid
     1775// because the spectral axis is the last one (so we can just
     1776// copy the first nAXesSub positions out)
     1777
     1778         pos = itData.pos().getFirst(nAxesSub);
     1779         sumSq(pos) += fac;
     1780//
     1781         itData.next();
     1782         itMask.next();
     1783      }
     1784   } else if (wtType==TSYS) {
     1785
     1786// We are going to average the data, weighted by 1/Tsys**2 for each pol, beam and IF.
     1787// So therefore we need to iterate through by spectrum (axis 3).  Although
     1788// Tsys is stored as a vector of length nChan, the values are replicated.
     1789// We will take a short cut and just use the value from the first channel
     1790// for now.
     1791//
     1792      VectorIterator<Float> itData(valuesIn, axis);
     1793      ReadOnlyVectorIterator<Float> itTSys(tSys, axis);
     1794      Float fac = 1.0;
     1795      IPosition pos(nAxesSub,0); 
     1796//
     1797      while (!itData.pastEnd()) {
     1798         Float t = itTSys.vector()[0];
     1799         fac = 1.0/t/t;
     1800
     1801// Scale data
     1802
     1803         itData.vector() *= fac;     // Writes back into 'dataIn'
     1804//
     1805// Accumulate Tsys  per if/pol/beam averaged over spectrum
     1806// This method to get pos2 from itData.pos() is only valid
     1807// because the spectral axis is the last one (so we can just
     1808// copy the first nAXesSub positions out)
     1809
     1810         pos = itData.pos().getFirst(nAxesSub);
     1811         tSysSqSum(pos) += fac;
     1812//
     1813         itData.next();
     1814         itTSys.next();
     1815      }
     1816   }
     1817
     1818// Accumulate sum of (possibly scaled) data
     1819
     1820   sum += dataIn;
     1821
     1822// Accumulate Tsys, time, and interval
     1823
     1824   tSysSum += tSys;
     1825   timeSum += time;
     1826   intSum += interval;
     1827   nAccum += 1;
     1828}
     1829
     1830
    17021831void SDMath::normalize(MaskedArray<Float>& sum,
    1703                         const Array<Float>& sumSq,
    1704                         const Array<Float>& nPts,
    1705                         WeightType wtType, Int axis,
    1706                         Int nAxesSub) const
     1832                       const Array<Float>& sumSq,
     1833                       const Array<Float>& tSysSqSum,
     1834                       const Array<Float>& nPts,
     1835                       Double intSum,
     1836                       WeightType wtType, Int axis,
     1837                       Int nAxesSub) const
    17071838{
    17081839   IPosition pos2(nAxesSub,0);
     
    17161847      MaskedArray<Float> t(nPts, (nPts>Float(0.0)));
    17171848      sum /= t;
     1849   } else if (wtType==TINT) {
     1850
     1851// Average by sum of Tint
     1852
     1853      sum /= Float(intSum);
    17181854   } else if (wtType==VAR) {
    17191855
     
    17291865      }
    17301866   } else if (wtType==TSYS) {
    1731    }
    1732 }
    1733 
    1734 
    1735 void SDMath::accumulate(Double& timeSum, Double& intSum, Int& nAccum,
    1736                         MaskedArray<Float>& sum, Array<Float>& sumSq,
    1737                         Array<Float>& nPts, Array<Float>& tSysSum,
    1738                         const Array<Float>& tSys, const Array<Float>& nInc,
    1739                         const Vector<Bool>& mask, Double time, Double interval,
    1740                         const Block<CountedPtr<SDMemTable> >& in,
    1741                         uInt iTab, uInt iRow, uInt axis,
    1742                         uInt nAxesSub, Bool useMask,
    1743                         WeightType wtType) const
    1744 {
    1745 
    1746 // Get data
    1747 
    1748    MaskedArray<Float> dataIn(in[iTab]->rowAsMaskedArray(iRow));
    1749    Array<Float>& valuesIn = dataIn.getRWArray();           // writable reference
    1750    const Array<Bool>& maskIn = dataIn.getMask();          // RO reference
    1751 //
    1752    if (wtType==NONE) {
    1753       const MaskedArray<Float> n(nInc,dataIn.getMask());
    1754       nPts += n;                               // Only accumulates where mask==T
    1755    } else if (wtType==VAR) {
    1756 
    1757 // We are going to average the data, weighted by the noise for each pol, beam and IF.
    1758 // So therefore we need to iterate through by spectrum (axis 3)
    1759 
    1760       VectorIterator<Float> itData(valuesIn, axis);
    1761       ReadOnlyVectorIterator<Bool> itMask(maskIn, axis);
    1762       Float fac = 1.0;
    1763       IPosition pos(nAxesSub,0); 
    1764 //
     1867   
     1868// Normalize each spectrum by sum(1/Tsys**2) where the pseudo
     1869// replication over channel for Tsys has been dropped.
     1870
     1871      Array<Float>& data = sum.getRWArray();
     1872      VectorIterator<Float> itData(data, axis);
    17651873      while (!itData.pastEnd()) {
    1766 
    1767 // Make MaskedArray of Vector, optionally apply OTF mask, and find scaling factor
    1768 
    1769         if (useMask) {
    1770            MaskedArray<Float> tmp(itData.vector(),mask&&itMask.vector());
    1771            fac = 1.0/variance(tmp);
    1772         } else {
    1773            MaskedArray<Float> tmp(itData.vector(),itMask.vector());
    1774            fac = 1.0/variance(tmp);
    1775         }
    1776 
    1777 // Scale data
    1778 
    1779         itData.vector() *= fac;     // Writes back into 'dataIn'
    1780 //
    1781 // Accumulate variance per if/pol/beam averaged over spectrum
    1782 // This method to get pos2 from itData.pos() is only valid
    1783 // because the spectral axis is the last one (so we can just
    1784 // copy the first nAXesSub positions out)
    1785 
    1786         pos = itData.pos().getFirst(nAxesSub);
    1787         sumSq(pos) += fac;
    1788 //
    1789         itData.next();
    1790         itMask.next();
     1874         pos2 = itData.pos().getFirst(nAxesSub);
     1875         itData.vector() /= tSysSqSum(pos2);
     1876         itData.next();
    17911877      }
    1792    } else if (wtType==TSYS) {
    1793    }
    1794 
    1795 // Accumulate sum of (possibly scaled) data
    1796 
    1797    sum += dataIn;
    1798 
    1799 // Accumulate Tsys, time, and interval
    1800 
    1801    tSysSum += tSys;
    1802    timeSum += time;
    1803    intSum += interval;
    1804    nAccum += 1;
    1805 }
     1878   }
     1879}
     1880
     1881
    18061882
    18071883
     
    18351911
    18361912
    1837 void SDMath::convertWeightString(WeightType& wtType, const String& weightStr) const
     1913void SDMath::convertWeightString(WeightType& wtType, const String& weightStr,
     1914                                 Bool listType) const
    18381915{
    18391916  String tStr(weightStr);
    18401917  tStr.upcase();
     1918  String msg;
    18411919  if (tStr.contains(String("NONE"))) {
    18421920     wtType = NONE;
     1921     msg = String("Weighting type selected : None");
    18431922  } else if (tStr.contains(String("VAR"))) {
    18441923     wtType = VAR;
     1924     msg = String("Weighting type selected : Variance");
     1925  } else if (tStr.contains(String("TINT"))) {
     1926     wtType = TINT;
     1927     msg = String("Weighting type selected : Tnt");
    18451928  } else if (tStr.contains(String("TSYS"))) {
    18461929     wtType = TSYS;
    1847      throw(AipsError("T_sys weighting not yet implemented"));
     1930     msg = String("Weighting type selected : Tsys");
    18481931  } else {
    1849     throw(AipsError("Unrecognized weighting type"));
    1850   }
     1932     msg = String("Weighting type selected : None");
     1933     throw(AipsError("Unrecognized weighting type"));
     1934  }
     1935//
     1936  if (listType) cout << msg << endl;
    18511937}
    18521938
  • trunk/src/SDMath.h

    r502 r518  
    136136// Weighting type for time averaging
    137137
    138   enum WeightType {NONE,VAR,TSYS};
     138  enum WeightType {NONE,VAR,TSYS,TINT};
    139139
    140140// Function to use accumulate data during time averaging
     
    143143                   casa::MaskedArray<casa::Float>& sum, casa::Array<casa::Float>& sumSq,
    144144                   casa::Array<casa::Float>& nPts, casa::Array<casa::Float>& tSysSum,
     145                   casa::Array<casa::Float>& tSysSqSum,
    145146                   const casa::Array<casa::Float>& tSys,  const casa::Array<casa::Float>& nInc,
    146147                   const casa::Vector<casa::Bool>& mask, casa::Double time, casa::Double interval,
     
    155156// Convert weight string to enum value
    156157
    157    void convertWeightString (WeightType& wt, const casa::String& weightStr) const;
     158   void convertWeightString (WeightType& wt, const casa::String& weightStr,
     159                             casa::Bool listType) const;
    158160
    159161// Convert interpolation type string
     
    233235  void normalize (casa::MaskedArray<casa::Float>& data,
    234236                  const casa::Array<casa::Float>& sumSq,
     237                  const casa::Array<casa::Float>& tSysSumSq,
    235238                  const casa::Array<casa::Float>& nPts,
     239                  casa::Double intSum,
    236240                  WeightType wtType, casa::Int axis, casa::Int nAxes) const;
    237241
Note: See TracChangeset for help on using the changeset viewer.