Changeset 130 for trunk


Ignore:
Timestamp:
12/20/04 15:44:42 (20 years ago)
Author:
kil064
Message:

Rewrite pretty much all of it to

  • use iterators rather than indexed for loops
  • loop over nrows

Fixed some minor bugs here and there
Added function 'statistics' to provide more statistics beyond standard deviation

Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/SDMath.cc

    r125 r130  
    3535#include <casa/Arrays/IPosition.h>
    3636#include <casa/Arrays/Array.h>
    37 #include <casa/Arrays/ArrayAccessor.h>
    38 #include <casa/Arrays/Slice.h>
     37#include <casa/Arrays/ArrayIter.h>
     38#include <casa/Arrays/VectorIter.h>
    3939#include <casa/Arrays/ArrayMath.h>
    4040#include <casa/Arrays/ArrayLogical.h>
    … …  
    4242#include <casa/Arrays/MaskArrMath.h>
    4343#include <casa/Arrays/MaskArrLogi.h>
    44 #include <casa/Arrays/VectorIter.h>
     44#include <casa/Exceptions.h>
    4545
    4646#include <tables/Tables/Table.h>
    … …  
    4848#include <tables/Tables/ArrayColumn.h>
    4949
    50 #include <scimath/Fitting.h>
    51 #include <scimath/Fitting/LinearFit.h>
    52 #include <scimath/Functionals/CompiledFunction.h>
    53 #include <images/Images/ImageUtilities.h>
     50#include <lattices/Lattices/LatticeUtilities.h>
     51#include <lattices/Lattices/RebinLattice.h>
    5452#include <coordinates/Coordinates/SpectralCoordinate.h>
    55 #include <scimath/Mathematics/AutoDiff.h>
    56 #include <scimath/Mathematics/AutoDiffMath.h>
     53#include <coordinates/Coordinates/CoordinateSystem.h>
     54#include <coordinates/Coordinates/CoordinateUtil.h>
    5755
    5856#include "MathUtils.h"
    … …  
    6664//using namespace asap::SDMath;
    6765
    68 CountedPtr<SDMemTable> SDMath::average(const CountedPtr<SDMemTable>& in) {
     66CountedPtr<SDMemTable> SDMath::average(const CountedPtr<SDMemTable>& in)
     67//
     68// Average all rows in Table in time
     69//
     70{
    6971  Table t = in->table();
    7072  ROArrayColumn<Float> tsys(t, "TSYS");
    … …  
    8587  Double inttime = 0.0;
    8688
     89// Loop over rows
     90
    8791  for (uInt i=0; i < t.nrow(); i++) {
    88     // data stuff
     92
     93// Get data and accumulate sums
     94
    8995    MaskedArray<Float> marr(in->rowAsMaskedArray(i));
    9096    outarr += marr;
    9197    MaskedArray<Float> n(narrinc,marr.getMask());
    9298    narr += n;
    93     // get
     99
     100// Accumulkate Tsys
     101
    94102    tsys.get(i, tsarr);// this is probably unneccessary as tsys should
    95103    outtsarr += tsarr; // be constant
    … …  
    100108    inttime += tmp;
    101109  }
    102   // averaging using mask
     110
     111// Average
     112
    103113  MaskedArray<Float> nma(narr,(narr > Float(0)));
    104114  outarr /= nma;
     115
     116// Create container and put
     117
    105118  Array<Bool> outflagsb = !(nma.getMask());
    106119  Array<uChar> outflags(outflagsb.shape());
    … …  
    127140CountedPtr<SDMemTable>
    128141SDMath::quotient(const CountedPtr<SDMemTable>& on,
    129                  const CountedPtr<SDMemTable>& off) {
     142                 const CountedPtr<SDMemTable>& off)
     143//
     144// Compute quotient spectrum
     145//
     146{
     147  const uInt nRows = on->nRow();
     148  if (off->nRow() != nRows) {
     149     throw (AipsError("Input Scan Tables must have the same number of rows"));
     150  }
     151
     152// Input Tables and columns
    130153
    131154  Table ton = on->table();
    … …  
    137160  ROArrayColumn<uInt> freqidc(ton, "FREQID");
    138161
    139   MaskedArray<Float> mon(on->rowAsMaskedArray(0));
    140   MaskedArray<Float> moff(off->rowAsMaskedArray(0));
    141   IPosition ipon = mon.shape();
    142   IPosition ipoff = moff.shape();
    143   Array<Float> tsarr;//(tsys.shape(0));
    144   tsys.get(0, tsarr);
    145   if (ipon != ipoff && ipon != tsarr.shape())
    146     cerr << "on/off not conformant" << endl;
    147 
    148   MaskedArray<Float> tmp = (mon-moff);
    149   Array<Float> out(tmp.getArray());
    150   out /= moff;
    151   out *= tsarr;
    152   Array<Bool> outflagsb = !(mon.getMask() && moff.getMask());
    153   Array<uChar> outflags(outflagsb.shape());
    154   convertArray(outflags,outflagsb);
    155   SDContainer sc = on->getSDContainer();
    156   sc.putTsys(tsarr);
    157   sc.scanid = 0;
    158   sc.putSpectrum(out);
    159   sc.putFlags(outflags);
     162// Output Table cloned from input
     163
    160164  SDMemTable* sdmt = new SDMemTable(*on, True);
    161   sdmt->putSDContainer(sc);
     165
     166// Loop over rows
     167
     168  for (uInt i=0; i<nRows; i++) {
     169     MaskedArray<Float> mon(on->rowAsMaskedArray(i));
     170     MaskedArray<Float> moff(off->rowAsMaskedArray(i));
     171     IPosition ipon = mon.shape();
     172     IPosition ipoff = moff.shape();
     173//
     174     Array<Float> tsarr; 
     175     tsys.get(i, tsarr);
     176     if (ipon != ipoff && ipon != tsarr.shape()) {
     177       throw(AipsError("on/off not conformant"));
     178     }
     179
     180// Compute quotient
     181
     182     MaskedArray<Float> tmp = (mon-moff);
     183     Array<Float> out(tmp.getArray());
     184     out /= moff;
     185     out *= tsarr;
     186     Array<Bool> outflagsb = !(mon.getMask() && moff.getMask());
     187     Array<uChar> outflags(outflagsb.shape());
     188     convertArray(outflags,outflagsb);
     189
     190// Fill container for this row
     191
     192     SDContainer sc = on->getSDContainer();
     193     sc.putTsys(tsarr);
     194     sc.scanid = 0;
     195     sc.putSpectrum(out);
     196     sc.putFlags(outflags);
     197
     198// Put new row in output Table
     199
     200     sdmt->putSDContainer(sc);
     201  }
     202//
    162203  return CountedPtr<SDMemTable>(sdmt);
    163204}
    164205
    165206CountedPtr<SDMemTable>
    166 SDMath::multiply(const CountedPtr<SDMemTable>& in, Float factor) {
     207SDMath::multiply(const CountedPtr<SDMemTable>& in, Float factor)
     208//
     209// Multiply values by factor
     210//
     211{
    167212  SDMemTable* sdmt = new SDMemTable(*in);
    168213  Table t = sdmt->table();
    … …  
    170215
    171216  for (uInt i=0; i < t.nrow(); i++) {
    172     // data stuff
    173217    MaskedArray<Float> marr(sdmt->rowAsMaskedArray(i));
    174218    marr *= factor;
    … …  
    179223
    180224CountedPtr<SDMemTable>
    181 SDMath::add(const CountedPtr<SDMemTable>& in, Float offset) {
     225SDMath::add(const CountedPtr<SDMemTable>& in, Float offset)
     226//
     227// Add offset to values
     228//
     229{
    182230  SDMemTable* sdmt = new SDMemTable(*in);
     231
    183232  Table t = sdmt->table();
    184233  ArrayColumn<Float> spec(t,"SPECTRA");
    185234
    186235  for (uInt i=0; i < t.nrow(); i++) {
    187     // data stuff
    188236    MaskedArray<Float> marr(sdmt->rowAsMaskedArray(i));
    189237    marr += offset;
    … …  
    194242
    195243
    196 // Start NEBK
    197 // SHow how to do above with VectorIterator
    198 
    199 //   uInt axis = 3;
    200 //   VectorIterator<Float> itData(arr, axis);
    201 //   ReadOnlyVectorIterator<Bool> itMask(barr, axis);
    202 //   Vector<Float> outv;
    203 //   while (!itData.pastEnd()) {
    204 //     SDMath::fit(outv, itData.vector(), itMask.vector()&&inmask, fitexpr);   // Should have function
    205 //     itData.vector() = outv;                                                 // to do in situ
    206 // //
    207 //     SDMath::fit(itData.vector(), itData.vector(), itMask.vector()&&inmask, fitexpr);   // Or maybe this works
    208 // //
    209 //     itData.next();
    210 //     itMask.next();
    211 //   }
    212 
    213 // End NEBK
    214 
    215244CountedPtr<SDMemTable>
    216 SDMath::hanning(const CountedPtr<SDMemTable>& in) {
    217  
    218   //IPosition ip = in->rowAsMaskedArray(0).shape();
     245SDMath::hanning(const CountedPtr<SDMemTable>& in)
     246//
     247// Hanning smooth each row
     248// Should Tsys be smoothed ?
     249//
     250{
    219251  SDMemTable* sdmt = new SDMemTable(*in,True);
     252
     253// Loop over rows in Table
     254
    220255  for (uInt ri=0; ri < in->nRow(); ++ri) {
     256
     257// Get data
    221258   
    222     MaskedArray<Float> marr(in->rowAsMaskedArray(ri));
    223 
     259    const MaskedArray<Float>& marr(in->rowAsMaskedArray(ri));
    224260    Array<Float> arr = marr.getArray();
    225261    Array<Bool> barr = marr.getMask();
    226     for (uInt i=0; i<in->nBeam();++i) {
    227         for (uInt j=0; j<in->nIF();++j) {
    228           for (uInt k=0; k<in->nPol();++k) {
    229             IPosition start(4,i,j,k,0);
    230             IPosition end(4,i,j,k,in->nChan()-1);
    231             Array<Float> subArr(arr(start,end));
    232             Array<Bool> subMask(barr(start,end));
    233             Vector<Float> outv;
    234             Vector<Bool> outm;
    235             Vector<Float> v(subArr.nonDegenerate());
    236             Vector<Bool> m(subMask.nonDegenerate());
    237             mathutil::hanning(outv,outm,v,m);
    238             ArrayAccessor<Float, Axis<0> > aa0(outv);
    239             ArrayAccessor<Bool, Axis<0> > ba0(outm);
    240             ArrayAccessor<Bool, Axis<3> > ba(subMask);
    241             for (ArrayAccessor<Float, Axis<3> > aa(subArr);
    242                  aa != aa.end();++aa) {
    243               (*aa) = (*aa0);
    244               (*ba) = (*ba0);
    245               aa0++;
    246               ba0++;
    247               ba++;
    248             }
    249           }
    250         }
     262
     263// Smooth along the channels axis
     264
     265    uInt axis = 3;
     266    VectorIterator<Float> itData(arr, axis);
     267    VectorIterator<Bool> itMask(barr, axis);
     268    Vector<Float> outv;
     269    Vector<Bool> outm;
     270    while (!itData.pastEnd()) {
     271       mathutil::hanning(outv, outm, itData.vector(), itMask.vector());
     272       itData.vector() = outv;                                                 
     273       itMask.vector() = outm;
     274//
     275       itData.next();
     276       itMask.next();
    251277    }
    252     //
    253    
    254 //     uInt axis = 3;
    255 //     VectorIterator<Float> itData(arr, axis);
    256 //     VectorIterator<Bool> itMask(barr, axis);
    257 //     Vector<Float> outv;
    258 //     Vector<Bool> outm;
    259 //     while (!itData.pastEnd()) {
    260 //       ::hanning(outv, outm, itData.vector(), itMask.vector());
    261 //       itData.vector() = outv;                                                 
    262 //       itMask.vector() = outm;
    263 //       //
    264 //       itData.next();
    265 //       itMask.next();
    266 //    }
    267    
    268     // End NEBK
    269    
    270     //
     278
     279// Create and put back
     280
    271281    Array<uChar> outflags(barr.shape());
    272282    convertArray(outflags,!barr);
    … …  
    279289}
    280290
    281 CountedPtr<SDMemTable>
    282 SDMath::averages(const Block<CountedPtr<SDMemTable> >& in,
    283                  const Vector<Bool>& mask) {
    284   IPosition ip = in[0]->rowAsMaskedArray(0).shape();
    285   Array<Float> arr(ip);
    286   Array<Bool> barr(ip);
    287   Double inttime = 0.0;
    288 
    289   uInt n = in[0]->nChan();
    290   for (uInt i=0; i<in[0]->nBeam();++i) {
    291     for (uInt j=0; j<in[0]->nIF();++j) {
    292       for (uInt k=0; k<in[0]->nPol();++k) {
    293         Float stdevsqsum = 0.0;
    294         Vector<Float> initvec(n);initvec = 0.0;
    295         Vector<Bool> initmask(n);initmask = True;
    296         MaskedArray<Float> outmarr(initvec,initmask);
    297         inttime = 0.0;
    298         for (uInt bi=0; bi< in.nelements(); ++bi) {
    299           MaskedArray<Float> marr(in[bi]->rowAsMaskedArray(0));
    300           inttime += in[bi]->getInterval();
    301           Array<Float> arr = marr.getArray();
    302           Array<Bool> barr = marr.getMask();
    303           IPosition start(4,i,j,k,0);
    304           IPosition end(4,i,j,k,n-1);
    305           Array<Float> subArr(arr(start,end));
    306           Array<Bool> subMask(barr(start,end));
    307           Vector<Float> outv;
    308           Vector<Bool> outm;
    309           Vector<Float> v(subArr.nonDegenerate());
    310           Vector<Bool> m(subMask.nonDegenerate());
    311           MaskedArray<Float> tmparr(v,m);
    312           MaskedArray<Float> tmparr2(tmparr(mask));
    313           Float stdvsq = pow(stddev(tmparr2),2);
    314           stdevsqsum+=1.0/stdvsq;
    315           tmparr /= stdvsq;
    316           outmarr += tmparr;
     291
     292
     293CountedPtr<SDMemTable> SDMath::averages(const Block<CountedPtr<SDMemTable> >& in,
     294                                        const Vector<Bool>& mask)
     295//
     296// Noise weighted averaging of spectra from many Tables.  Tables can have different
     297// number of rows. 
     298//
     299{
     300
     301// Setup
     302
     303  const uInt axis = 3;                                 // Spectral axis
     304  IPosition shp = in[0]->rowAsMaskedArray(0).shape();
     305  Array<Float> arr(shp);
     306  Array<Bool> barr(shp);
     307  Double sumInterval = 0.0;
     308  const Bool useMask = (mask.nelements() == shp(axis));
     309
     310// Create data accumulation MaskedArray. We accumulate for each
     311// channel,if,pol,beam
     312
     313  Array<Float> zero(shp); zero=0.0;
     314  Array<Bool> good(shp); good = True;
     315  MaskedArray<Float> sum(zero,good);
     316
     317// Create accumulation Array for variance. We accumulate for
     318// each if,pol,beam, but average over channel
     319
     320  const uInt nAxesSub = shp.nelements() - 1;
     321  IPosition shp2(nAxesSub);
     322  for (uInt i=0,j=0; i<(nAxesSub+1); i++) {
     323     if (i!=axis) {
     324       shp2(j) = shp(i);
     325       j++;
     326     }
     327  }
     328  Array<Float> sumSq(shp2);
     329  sumSq = 0.0;
     330  IPosition pos2(nAxesSub,0);                        // FOr indexing
     331// 
     332  Float fac = 1.0;
     333  const uInt nTables = in.nelements();
     334  for (uInt iTab=0; iTab<nTables; iTab++) {
     335     const uInt nRows = in[iTab]->nRow();
     336     sumInterval += nRows * in[iTab]->getInterval();   // Sum of time intervals
     337//
     338     for (uInt iRow=0; iRow<nRows; iRow++) {
     339
     340// Check conforms
     341
     342        IPosition shp2 = in[iTab]->rowAsMaskedArray(iRow).shape();
     343        if (!shp.isEqual(shp2)) {
     344           throw (AipsError("Shapes for all rows must be the same"));
    317345        }
    318         outmarr /= stdevsqsum;
    319         Array<Float> tarr(outmarr.getArray());
    320         Array<Bool> tbarr(outmarr.getMask());
    321         IPosition start(4,i,j,k,0);
    322         IPosition end(4,i,j,k,n-1);
    323         Array<Float> subArr(arr(start,end));
    324         Array<Bool> subMask(barr(start,end));
    325         ArrayAccessor<Float, Axis<0> > aa0(tarr);
    326         ArrayAccessor<Bool, Axis<0> > ba0(tbarr);
    327         ArrayAccessor<Bool, Axis<3> > ba(subMask);
    328         for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
    329           (*aa) = (*aa0);
    330           (*ba) = (*ba0);
    331           aa0++;
    332           ba0++;
    333           ba++;
    334         }
    335       }
     346
     347// Get data and deconstruct
     348 
     349        MaskedArray<Float> marr(in[iTab]->rowAsMaskedArray(iRow));
     350        Array<Float>& arr = marr.getRWArray();                     // writable reference
     351        const Array<Bool>& barr = marr.getMask();                  // RO reference
     352
     353// We are going to average the data, weighted by the noise for each
     354// pol, beam and IF. So therefore we need to iterate through by
     355// spectra (axis 3)
     356
     357        VectorIterator<Float> itData(arr, axis);
     358        ReadOnlyVectorIterator<Bool> itMask(barr, axis);
     359        while (!itData.pastEnd()) {
     360
     361// Make MaskedArray of Vector, optionally apply OTF mask, and find scaling factor
     362
     363          if (useMask) {
     364             MaskedArray<Float> tmp(itData.vector(),mask&&itMask.vector());
     365             fac = 1.0/variance(tmp);
     366          } else {
     367             MaskedArray<Float> tmp(itData.vector(),itMask.vector());
     368             fac = 1.0/variance(tmp);
     369          }
     370
     371// Scale data
     372
     373          itData.vector() *= fac;
     374
     375// Accumulate variance per if/pol/beam averaged over spectrum
     376// This method to get pos2 from itData.pos() is only valid
     377// because the spectral axis is the last one (so we can just
     378// copy the first nAXesSub positions out)
     379
     380          pos2 = itData.pos().getFirst(nAxesSub);
     381          sumSq(pos2) += fac;
     382//
     383          itData.next();
     384          itMask.next();
     385        }   
     386
     387// Accumulate sums
     388
     389       sum += marr;
    336390    }
    337391  }
    338   Array<uChar> outflags(barr.shape());
    339   convertArray(outflags,!barr);
    340   SDContainer sc = in[0]->getSDContainer();
    341   sc.putSpectrum(arr);
     392
     393// Normalize by the sum of the 1/var. 
     394
     395  Array<Float>& data = sum.getRWArray();   
     396  VectorIterator<Float> itData(data, axis);
     397  while (!itData.pastEnd()) {
     398     pos2 = itData.pos().getFirst(nAxesSub);           // See comments above
     399     itData.vector() /= sumSq(pos2);
     400     itData.next();
     401  }   
     402
     403// Create and fill output
     404
     405  Array<uChar> outflags(shp);
     406  convertArray(outflags,!(sum.getMask()));
     407//
     408  SDContainer sc = in[0]->getSDContainer();     // CLone from first container of first Table
     409  sc.putSpectrum(data);
    342410  sc.putFlags(outflags);
    343   sc.interval = inttime;
    344   SDMemTable* sdmt = new SDMemTable(*in[0],True);
     411  sc.interval = sumInterval;
     412//
     413  SDMemTable* sdmt = new SDMemTable(*in[0],True);  // CLone from first Table
    345414  sdmt->putSDContainer(sc);
    346415  return CountedPtr<SDMemTable>(sdmt);
    347416}
     417
    348418
    349419CountedPtr<SDMemTable>
    350420SDMath::averagePol(const CountedPtr<SDMemTable>& in,
    351                    const Vector<Bool>& mask) {
    352   MaskedArray<Float> marr(in->rowAsMaskedArray(0));
    353   uInt n = in->nChan();
    354   IPosition ip = marr.shape();
    355   Array<Float> arr = marr.getArray();
    356   Array<Bool> barr = marr.getMask();
    357   for (uInt i=0; i<in->nBeam();++i) {
    358     for (uInt j=0; j<in->nIF();++j) {
    359       Float stdevsqsum = 0.0;
    360       Vector<Float> initvec(n);initvec = 0.0;
    361       Vector<Bool> initmask(n);initmask = True;
    362       MaskedArray<Float> outmarr(initvec,initmask);
    363       for (uInt k=0; k<in->nPol();++k) {
    364         IPosition start(4,i,j,k,0);
    365         IPosition end(4,i,j,k,in->nChan()-1);
    366         Array<Float> subArr(arr(start,end));
    367         Array<Bool> subMask(barr(start,end));
    368         Vector<Float> outv;
    369         Vector<Bool> outm;
    370         Vector<Float> v(subArr.nonDegenerate());
    371         Vector<Bool> m(subMask.nonDegenerate());
    372         MaskedArray<Float> tmparr(v,m);
    373         MaskedArray<Float> tmparr2(tmparr(mask));
    374         Float stdvsq = pow(stddev(tmparr2),2);
    375         stdevsqsum+=1.0/stdvsq;
    376         tmparr /= stdvsq;
    377         outmarr += tmparr;
    378       }
    379       outmarr /= stdevsqsum;
    380       Array<Float> tarr(outmarr.getArray());
    381       Array<Bool> tbarr(outmarr.getMask());
    382       // write averaged pol into all pols - fix up to reform array
    383       for (uInt k=0; k<in->nPol();++k) {
    384         IPosition start(4,i,j,k,0);
    385         IPosition end(4,i,j,k,n-1);
    386         Array<Float> subArr(arr(start,end));
    387         Array<Bool> subMask(barr(start,end));
    388         ArrayAccessor<Float, Axis<0> > aa0(tarr);
    389         ArrayAccessor<Bool, Axis<0> > ba0(tbarr);
    390         ArrayAccessor<Bool, Axis<3> > ba(subMask);
    391         for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
    392           (*aa) = (*aa0);
    393           (*ba) = (*ba0);
    394           aa0++;
    395           ba0++;
    396           ba++;
     421                   const Vector<Bool>& mask)
     422{
     423   const uInt nRows = in->nRow();
     424   const uInt axis = 3;                        // Spectrum
     425   const IPosition axes(2, 2, 3);              // pol-channel plane
     426
     427// Create output Table
     428
     429  SDMemTable* sdmt = new SDMemTable(*in, True);
     430
     431// Loop over rows
     432
     433   for (uInt iRow=0; iRow<nRows; iRow++) {
     434
     435// Get data for this row
     436
     437      MaskedArray<Float> marr(in->rowAsMaskedArray(iRow));
     438      Array<Float>& arr = marr.getRWArray();
     439      const Array<Bool>& barr = marr.getMask();
     440//
     441      IPosition shp = marr.shape();
     442      const Bool useMask = (mask.nelements() == shp(axis));
     443      const uInt nChan = shp(axis);
     444
     445// Make iterators to iterate by pol-channel planes
     446
     447     ArrayIterator<Float> itDataPlane(arr, axes);
     448     ReadOnlyArrayIterator<Bool> itMaskPlane(barr, axes);
     449
     450// Accumulations
     451
     452     Float fac = 0.0;
     453     Vector<Float> vecSum(nChan,0.0);
     454
     455// Iterate by plane
     456
     457     while (!itDataPlane.pastEnd()) {
     458
     459// Iterate through pol-channel plane by spectrum
     460
     461        Vector<Float> t1(nChan); t1 = 0.0;
     462        Vector<Bool> t2(nChan); t2 = True;
     463        MaskedArray<Float> vecSum(t1,t2);
     464        Float varSum = 0.0;
     465        {
     466           ReadOnlyVectorIterator<Float> itDataVec(itDataPlane.array(), 1);
     467           ReadOnlyVectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
     468           while (!itDataVec.pastEnd()) {     
     469
     470// Create MA of data & mask (optionally including OTF mask) and  get variance
     471
     472              if (useMask) {
     473                 const MaskedArray<Float> spec(itDataVec.vector(),mask&&itMaskVec.vector());
     474                 fac = 1.0 / variance(spec);
     475              } else {
     476                 const MaskedArray<Float> spec(itDataVec.vector(),itMaskVec.vector());
     477                 fac = 1.0 / variance(spec);
     478              }
     479
     480// Normalize spectrum (without OTF mask) and accumulate
     481
     482              const MaskedArray<Float> spec(fac*itDataVec.vector(), itMaskVec.vector());
     483              vecSum += spec;
     484              varSum += fac;
     485
     486// Next
     487
     488              itDataVec.next();
     489              itMaskVec.next();
     490           }
    397491        }
    398       }
    399     }
    400   }
    401 
    402   Array<uChar> outflags(barr.shape());
    403   convertArray(outflags,!barr);
    404   SDContainer sc = in->getSDContainer();
    405   sc.putSpectrum(arr);
    406   sc.putFlags(outflags);
    407   SDMemTable* sdmt = new SDMemTable(*in,True);
    408   sdmt->putSDContainer(sc);
    409   return CountedPtr<SDMemTable>(sdmt);
    410 }
    411 
    412 
    413 Float SDMath::rms(const CountedPtr<SDMemTable>& in,
    414                          const std::vector<bool>& mask) {
    415   Float rmsval;
    416   Vector<Bool> msk(mask);
    417   IPosition ip = in->rowAsMaskedArray(0).shape();
    418   MaskedArray<Float> marr(in->rowAsMaskedArray(0));
    419 
    420   Array<Float> arr = marr.getArray();
    421   Array<Bool> barr = marr.getMask();
    422   uInt i,j,k;
    423   i = in->getBeam();
    424   j = in->getIF();
    425   k = in->getPol();
    426   IPosition start(4,i,j,k,0);
    427   IPosition end(4,i,j,k,in->nChan()-1);
    428   Array<Float> subArr(arr(start,end));
    429   Array<Bool> subMask(barr(start,end));
    430   Array<Float> v(subArr.nonDegenerate());
    431   Array<Bool> m(subMask.nonDegenerate());
    432   MaskedArray<Float> tmp;
    433   if (msk.nelements() == m.nelements() ) {
    434     tmp.setData(v,m&&msk);
    435   } else {
    436     tmp.setData(v,m);
    437   }
    438   rmsval = stddev(tmp);
    439   return rmsval;
    440 }
     492
     493// Normalize summed spectrum
     494
     495        vecSum /= varSum;
     496
     497// We have formed the weighted averaged spectrum from all polarizations
     498// for this beam and IF.  Now replicate the spectrum to all polarizations
     499
     500        {
     501           VectorIterator<Float> itDataVec(itDataPlane.array(), 1);  // Writes back into 'arr'
     502           const Vector<Float>& vecSumData = vecSum.getArray();      // It *is* a Vector
     503//
     504           while (!itDataVec.pastEnd()) {     
     505              itDataVec.vector() = vecSumData;
     506              itDataVec.next();
     507           }
     508        }
     509
     510// Step to next beam/IF combination
     511
     512        itDataPlane.next();
     513        itMaskPlane.next();
     514     }
     515
     516// Generate output container and write it to output table
     517
     518     SDContainer sc = in->getSDContainer();
     519     Array<uChar> outflags(barr.shape());
     520     convertArray(outflags,!barr);
     521     sc.putSpectrum(arr);
     522     sc.putFlags(outflags);
     523     sdmt->putSDContainer(sc);
     524   }
     525//
     526  return CountedPtr<SDMemTable>(sdmt);
     527}
     528
    441529
    442530CountedPtr<SDMemTable> SDMath::bin(const CountedPtr<SDMemTable>& in,
    443                                    Int width) {
    444 
     531                                   Int width)
     532{
    445533  SDHeader sh = in->getSDHeader();
    446534  SDMemTable* sdmt = new SDMemTable(*in,True);
    447535
    448   MaskedArray<Float> tmpmarr(in->rowAsMaskedArray(0));
    449   MaskedArray<Float> tmpmarr2;
    450   SpectralCoordinate outcoord;
    451   IPosition ip;
     536// Bin up SpectralCoordinates
     537
     538  IPosition factors(1);
     539  factors(0) = width;
    452540  for (uInt j=0; j<in->nCoordinates(); ++j) {
    453     SpectralCoordinate coord(in->getCoordinate(j));
    454     ImageUtilities::bin(tmpmarr2, outcoord, tmpmarr, coord, 3, width);
    455     ip = tmpmarr2.shape();
    456     sdmt->setCoordinate(outcoord,j);
    457   }
    458   sh.nchan = ip(3);
     541    CoordinateSystem cSys;
     542    cSys.addCoordinate(in->getCoordinate(j));
     543    CoordinateSystem cSysBin =
     544      CoordinateUtil::makeBinnedCoordinateSystem (factors, cSys, False);
     545//
     546    SpectralCoordinate sCBin = cSysBin.spectralCoordinate(0);
     547    sdmt->setCoordinate(sCBin, j);
     548  }
     549
     550// Use RebinLattice to find shape
     551
     552  IPosition shapeIn(1,sh.nchan);
     553  IPosition shapeOut = RebinLattice<Float>::rebinShape (shapeIn, factors);
     554  sh.nchan = shapeOut(0);
    459555  sdmt->putSDHeader(sh);
    460556
    461   SpectralCoordinate tmpcoord2,tmpcoord;
     557
     558// Loop over rows and bin along channel axis
     559 
     560  const uInt axis = 3;
    462561  for (uInt i=0; i < in->nRow(); ++i) {
     562    SDContainer sc = in->getSDContainer(i);
     563//
     564    Array<Float> tSys(sc.getTsys());                           // Get it out before sc changes shape
     565
     566// Bin up spectrum
     567
    463568    MaskedArray<Float> marr(in->rowAsMaskedArray(i));
    464     SDContainer sc = in->getSDContainer(i);
    465     Array<Float> arr = marr.getArray();
    466     Array<Bool> barr = marr.getMask();
    467569    MaskedArray<Float> marrout;
    468     ImageUtilities::bin(marrout, tmpcoord2, marr, tmpcoord, 3, width);
     570    LatticeUtilities::bin(marrout, marr, axis, width);
     571
     572// Put back the binned data and flags
     573
    469574    IPosition ip2 = marrout.shape();
    470575    sc.resize(ip2);
    471576    sc.putSpectrum(marrout.getArray());
     577//
    472578    Array<uChar> outflags(ip2);
    473579    convertArray(outflags,!(marrout.getMask()));
    474580    sc.putFlags(outflags);
    475     MaskedArray<Float> tsys,tsysout;
    476     Array<Bool> dummy(ip);dummy = True;
    477     tsys.setData(sc.getTsys(),dummy);
    478     ImageUtilities::bin(tsysout, tmpcoord2, marr, tmpcoord, 3, width);
    479     sc.putTsys(tsysout.getArray());
     581
     582// Bin up Tsys. 
     583
     584    Array<Bool> allGood(tSys.shape(),True);
     585    MaskedArray<Float> tSysIn(tSys, allGood, True);
     586//
     587    MaskedArray<Float> tSysOut;   
     588    LatticeUtilities::bin(tSysOut, tSysIn, axis, width);
     589    sc.putTsys(tSysOut.getArray());
    480590    sdmt->putSDContainer(sc);
    481591  }
    482592  return CountedPtr<SDMemTable>(sdmt);
    483593}
     594
     595
     596
     597std::vector<float> SDMath::statistic (const CountedPtr<SDMemTable>& in,
     598                                       const std::vector<bool>& mask,
     599                                       const std::string& which)
     600//
     601// Perhaps iteration over pol/beam/if should be in here
     602// and inside the nrow iteration ?
     603//
     604{
     605  const uInt nRow = in->nRow();
     606  std::vector<float> result(nRow);
     607  Vector<Bool> msk(mask);
     608
     609// Specify cursor location
     610
     611  uInt i = in->getBeam();
     612  uInt j = in->getIF();
     613  uInt k = in->getPol();
     614  IPosition start(4,i,j,k,0);
     615  IPosition end(4,i,j,k,in->nChan()-1);
     616
     617// Loop over rows
     618
     619  const uInt nEl = msk.nelements();
     620  for (uInt ii=0; ii < in->nRow(); ++ii) {
     621
     622// Get row and deconstruct
     623
     624     MaskedArray<Float> marr(in->rowAsMaskedArray(ii));
     625     Array<Float> arr = marr.getArray();
     626     Array<Bool> barr = marr.getMask();
     627
     628// Access desired piece of data
     629
     630     Array<Float> v((arr(start,end)).nonDegenerate());
     631     Array<Bool> m((barr(start,end)).nonDegenerate());
     632
     633// Apply OTF mask
     634
     635     MaskedArray<Float> tmp;
     636     if (m.nelements()==nEl) {
     637       tmp.setData(v,m&&msk);
     638     } else {
     639       tmp.setData(v,m);
     640     }
     641
     642// Get statistic
     643
     644     result[ii] = SDMath::theStatistic(which, tmp);
     645  }
     646//
     647  return result;
     648}
     649
     650
     651float SDMath::theStatistic(const std::string& which,  const casa::MaskedArray<Float>& data)
     652{
     653   String str(which);
     654   str.upcase();
     655   if (str.contains(String("MIN"))) {
     656      return min(data);
     657   } else if (str.contains(String("MAX"))) {
     658      return max(data);
     659   } else if (str.contains(String("SUMSQ"))) {
     660      return sumsquares(data);
     661   } else if (str.contains(String("SUM"))) {
     662      return sum(data);
     663   } else if (str.contains(String("MEAN"))) {
     664      return mean(data);
     665   } else if (str.contains(String("VAR"))) {
     666      return variance(data);
     667   } else if (str.contains(String("STDDEV"))) {
     668      return stddev(data);
     669   } else if (str.contains(String("AVDEV"))) {
     670      return avdev(data);
     671   } else if (str.contains(String("RMS"))) {
     672      uInt n = data.nelementsValid();
     673      return sqrt(sumsquares(data)/n);
     674   } else if (str.contains(String("MED"))) {
     675      return median(data);
     676   }
     677}

  • trunk/src/SDMath.h

    r125 r130  
    3434#include <string>
    3535#include <vector>
     36#include <casa/aips.h>
    3637#include <casa/Utilities/CountedPtr.h>
    3738
    … …  
    5960  averagePol(const casa::CountedPtr<SDMemTable>& in, const casa::Vector<casa::Bool>& mask);
    6061
    61   casa::Float rms(const casa::CountedPtr<SDMemTable>& in,
    62                    const std::vector<bool>& mask);
     62  std::vector<float> statistic(const casa::CountedPtr<SDMemTable>& in,
     63                                const std::vector<bool>& mask, const std::string& which);
    6364 
    6465  casa::CountedPtr<SDMemTable> bin(const casa::CountedPtr<SDMemTable>& in,
    6566                             casa::Int width);
     67
     68// private (not actually...)
     69
     70  float theStatistic(const std::string& which,  const casa::MaskedArray<casa::Float>& data);
     71 
    6672};
    6773
Note: See TracChangeset for help on using the changeset viewer.