Changeset 85 for trunk/src/SDMath.cc


Ignore:
Timestamp:
11/12/04 12:03:40 (20 years ago)
Author:
mar637
Message:

bin function now bins all rows in the table. ALso fixed a couple of bugs.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/SDMath.cc

    r83 r85  
    6666CountedPtr<SDMemTable> SDMath::average(const CountedPtr<SDMemTable>& in) {
    6767  Table t = in->table();
    68   ROArrayColumn<Float> tsys(t, "TSYS"); 
     68  ROArrayColumn<Float> tsys(t, "TSYS");
    6969  ROScalarColumn<Double> mjd(t, "TIME");
    7070  ROScalarColumn<String> srcn(t, "SRCNAME");
    … …  
    8989    MaskedArray<Float> n(narrinc,marr.getMask());
    9090    narr += n;
    91     // get 
     91    // get
    9292    tsys.get(i, tsarr);// this is probably unneccessary as tsys should
    9393    outtsarr += tsarr; // be constant
    … …  
    117117  sc.scanid = 0;
    118118  sc.putSpectrum(outarr);
    119   sc.putFlags(outflags); 
     119  sc.putFlags(outflags);
    120120  SDMemTable* sdmt = new SDMemTable(*in,True);
    121121  sdmt->putSDContainer(sc);
    … …  
    123123}
    124124
    125 CountedPtr<SDMemTable> 
    126 SDMath::quotient(const CountedPtr<SDMemTable>& on, 
    127                 const CountedPtr<SDMemTable>& off) {
    128  
     125CountedPtr<SDMemTable>
     126SDMath::quotient(const CountedPtr<SDMemTable>& on,
     127                const CountedPtr<SDMemTable>& off) {
     128
    129129  Table ton = on->table();
    130130  Table toff = off->table();
    131   ROArrayColumn<Float> tsys(toff, "TSYS"); 
     131  ROArrayColumn<Float> tsys(toff, "TSYS");
    132132  ROScalarColumn<Double> mjd(ton, "TIME");
    133133  ROScalarColumn<Double> integr(ton, "INTERVAL");
    … …  
    141141  Array<Float> tsarr;//(tsys.shape(0));
    142142  tsys.get(0, tsarr);
    143   if (ipon != ipoff && ipon != tsarr.shape()) 
     143  if (ipon != ipoff && ipon != tsarr.shape())
    144144    cerr << "on/off not conformant" << endl;
    145  
     145
    146146  MaskedArray<Float> tmp = (mon-moff);
    147   Array<Float> out(tmp.getArray()); 
     147  Array<Float> out(tmp.getArray());
    148148  out /= moff;
    149149  out *= tsarr;
    … …  
    161161}
    162162
    163 CountedPtr<SDMemTable> 
     163CountedPtr<SDMemTable>
    164164SDMath::multiply(const CountedPtr<SDMemTable>& in, Float factor) {
    165165  SDMemTable* sdmt = new SDMemTable(*in);
    … …  
    176176}
    177177
    178 bool SDMath::fit(Vector<Float>& thefit, const Vector<Float>& data, 
    179                 const Vector<Bool>& mask,
    180                 const std::string& fitexpr) {
     178bool SDMath::fit(Vector<Float>& thefit, const Vector<Float>& data,
     179                const Vector<Bool>& mask,
     180                const std::string& fitexpr) {
    181181
    182182  LinearFit<Float> fitter;
    … …  
    194194}
    195195
    196 CountedPtr<SDMemTable> 
    197 SDMath::baseline(const CountedPtr<SDMemTable>& in, 
    198                 const std::string& fitexpr,
    199                 const std::vector<bool>& mask) {
    200  
     196CountedPtr<SDMemTable>
     197SDMath::baseline(const CountedPtr<SDMemTable>& in,
     198                const std::string& fitexpr,
     199                const std::vector<bool>& mask) {
     200
    201201  IPosition ip = in->rowAsMaskedArray(0).shape();
    202202  SDContainer sc = in->getSDContainer();
    … …  
    212212    for (uInt j=0; j<in->nIF();++j) {
    213213      for (uInt k=0; k<in->nPol();++k) {
    214         IPosition start(4,i,j,k,0);
    215         IPosition end(4,i,j,k,in->nChan()-1);
    216         Array<Float> subArr(arr(start,end));
    217         Array<Bool> subMask(barr(start,end));
    218         Vector<Float> outv;
    219         Vector<Float> v(subArr.nonDegenerate());
    220         Vector<Bool> m(subMask.nonDegenerate());
    221         cout << "\t Polarisation " << k << "\t";
    222         SDMath::fit(outv, v, m&&inmask, fitexpr);
    223         ArrayAccessor<Float, Axis<0> > aa0(outv);
    224         for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
    225           (*aa) = (*aa0);
    226           aa0++;
    227         }
     214        IPosition start(4,i,j,k,0);
     215        IPosition end(4,i,j,k,in->nChan()-1);
     216        Array<Float> subArr(arr(start,end));
     217        Array<Bool> subMask(barr(start,end));
     218        Vector<Float> outv;
     219        Vector<Float> v(subArr.nonDegenerate());
     220        Vector<Bool> m(subMask.nonDegenerate());
     221        cout << "\t Polarisation " << k << "\t";
     222        SDMath::fit(outv, v, m&&inmask, fitexpr);
     223        ArrayAccessor<Float, Axis<0> > aa0(outv);
     224        for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
     225          (*aa) = (*aa0);
     226          aa0++;
     227        }
    228228      }
    229229    }
    … …  
    239239
    240240
    241 CountedPtr<SDMemTable> 
     241CountedPtr<SDMemTable>
    242242SDMath::hanning(const CountedPtr<SDMemTable>& in) {
    243243
    244   IPosition ip = in->rowAsMaskedArray(0).shape();
    245   MaskedArray<Float> marr(in->rowAsMaskedArray(0));
    246 
    247   Array<Float> arr = marr.getArray();
    248   Array<Bool> barr = marr.getMask();
    249   for (uInt i=0; i<in->nBeam();++i) {
    250     for (uInt j=0; j<in->nIF();++j) {
    251       for (uInt k=0; k<in->nPol();++k) {
    252         IPosition start(4,i,j,k,0);
    253         IPosition end(4,i,j,k,in->nChan()-1);
    254         Array<Float> subArr(arr(start,end));
    255         Array<Bool> subMask(barr(start,end));
    256         Vector<Float> outv;
    257         Vector<Bool> outm;
    258         Vector<Float> v(subArr.nonDegenerate());
    259         Vector<Bool> m(subMask.nonDegenerate());
    260         ::hanning(outv,outm,v,m);
    261         ArrayAccessor<Float, Axis<0> > aa0(outv);       
    262         ArrayAccessor<Bool, Axis<0> > ba0(outm);       
    263         ArrayAccessor<Bool, Axis<3> > ba(subMask);     
    264         for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
    265           (*aa) = (*aa0);
    266           (*ba) = (*ba0);
    267           aa0++;
    268           ba0++;
    269           ba++;
    270         }
    271       }
     244  //IPosition ip = in->rowAsMaskedArray(0).shape();
     245  SDMemTable* sdmt = new SDMemTable(*in,True);
     246  for (uInt ri=0; ri < in->nRow(); ++ri) {
     247
     248    MaskedArray<Float> marr(in->rowAsMaskedArray(ri));
     249
     250    Array<Float> arr = marr.getArray();
     251    Array<Bool> barr = marr.getMask();
     252    for (uInt i=0; i<in->nBeam();++i) {
     253        for (uInt j=0; j<in->nIF();++j) {
     254        for (uInt k=0; k<in->nPol();++k) {
     255            IPosition start(4,i,j,k,0);
     256            IPosition end(4,i,j,k,in->nChan()-1);
     257            Array<Float> subArr(arr(start,end));
     258            Array<Bool> subMask(barr(start,end));
     259            Vector<Float> outv;
     260            Vector<Bool> outm;
     261            Vector<Float> v(subArr.nonDegenerate());
     262            Vector<Bool> m(subMask.nonDegenerate());
     263            ::hanning(outv,outm,v,m);
     264            ArrayAccessor<Float, Axis<0> > aa0(outv);
     265            ArrayAccessor<Bool, Axis<0> > ba0(outm);
     266            ArrayAccessor<Bool, Axis<3> > ba(subMask);
     267            for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
     268            (*aa) = (*aa0);
     269            (*ba) = (*ba0);
     270            aa0++;
     271            ba0++;
     272            ba++;
     273            }
     274        }
     275        }
    272276    }
    273   }
    274   Array<uChar> outflags(barr.shape());
    275   convertArray(outflags,!barr);
    276   SDContainer sc = in->getSDContainer();
    277   sc.putSpectrum(arr);
    278   sc.putFlags(outflags);
    279   SDMemTable* sdmt = new SDMemTable(*in,True);
    280   sdmt->putSDContainer(sc);
    281   return CountedPtr<SDMemTable>(sdmt);
    282 }
    283 
    284 CountedPtr<SDMemTable>
     277    Array<uChar> outflags(barr.shape());
     278    convertArray(outflags,!barr);
     279    SDContainer sc = in->getSDContainer(ri);
     280    sc.putSpectrum(arr);
     281    sc.putFlags(outflags);
     282    sdmt->putSDContainer(sc);
     283  }
     284  return CountedPtr<SDMemTable>(sdmt);
     285}
     286
     287CountedPtr<SDMemTable>
    285288SDMath::averages(const Block<CountedPtr<SDMemTable> >& in,
    286                 const Vector<Bool>& mask) {
     289                const Vector<Bool>& mask) {
    287290  IPosition ip = in[0]->rowAsMaskedArray(0).shape();
    288291  Array<Float> arr(ip);
    289292  Array<Bool> barr(ip);
    290293  Double inttime = 0.0;
    291  
     294
    292295  uInt n = in[0]->nChan();
    293296  for (uInt i=0; i<in[0]->nBeam();++i) {
    294297    for (uInt j=0; j<in[0]->nIF();++j) {
    295298      for (uInt k=0; k<in[0]->nPol();++k) {
    296         Float stdevsqsum = 0.0;
    297         Vector<Float> initvec(n);initvec = 0.0;
    298         Vector<Bool> initmask(n);initmask = True;
    299         MaskedArray<Float> outmarr(initvec,initmask);
    300         for (uInt bi=0; bi< in.nelements(); ++bi) {
    301           MaskedArray<Float> marr(in[bi]->rowAsMaskedArray(0));
    302           inttime += in[bi]->getInterval();
    303           Array<Float> arr = marr.getArray();
    304           Array<Bool> barr = marr.getMask();
    305           IPosition start(4,i,j,k,0);
    306           IPosition end(4,i,j,k,n-1);
    307           Array<Float> subArr(arr(start,end));
    308           Array<Bool> subMask(barr(start,end));
    309           Vector<Float> outv;
    310           Vector<Bool> outm;
    311           Vector<Float> v(subArr.nonDegenerate());
    312           Vector<Bool> m(subMask.nonDegenerate());
    313           MaskedArray<Float> tmparr(v,m);
    314           MaskedArray<Float> tmparr2(tmparr(mask));
    315           Float stdvsq = pow(stddev(tmparr2),2);
    316           stdevsqsum+=1.0/stdvsq;
    317           tmparr /= stdvsq;
    318           outmarr += tmparr;
    319         }
    320         outmarr /= stdevsqsum;
    321         Array<Float> tarr(outmarr.getArray());
    322         Array<Bool> tbarr(outmarr.getMask());
    323         IPosition start(4,i,j,k,0);
    324         IPosition end(4,i,j,k,n-1);
    325         Array<Float> subArr(arr(start,end));
    326         Array<Bool> subMask(barr(start,end));
    327         ArrayAccessor<Float, Axis<0> > aa0(tarr);       
    328         ArrayAccessor<Bool, Axis<0> > ba0(tbarr);       
    329         ArrayAccessor<Bool, Axis<3> > ba(subMask);     
    330         for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
    331           (*aa) = (*aa0);
    332           (*ba) = (*ba0);
    333           aa0++;
    334           ba0++;
    335           ba++;
    336         }
     299        Float stdevsqsum = 0.0;
     300        Vector<Float> initvec(n);initvec = 0.0;
     301        Vector<Bool> initmask(n);initmask = True;
     302        MaskedArray<Float> outmarr(initvec,initmask);
     303        inttime = 0.0;
     304        for (uInt bi=0; bi< in.nelements(); ++bi) {
     305          MaskedArray<Float> marr(in[bi]->rowAsMaskedArray(0));
     306          inttime += in[bi]->getInterval();
     307          Array<Float> arr = marr.getArray();
     308          Array<Bool> barr = marr.getMask();
     309          IPosition start(4,i,j,k,0);
     310          IPosition end(4,i,j,k,n-1);
     311          Array<Float> subArr(arr(start,end));
     312          Array<Bool> subMask(barr(start,end));
     313          Vector<Float> outv;
     314          Vector<Bool> outm;
     315          Vector<Float> v(subArr.nonDegenerate());
     316          Vector<Bool> m(subMask.nonDegenerate());
     317          MaskedArray<Float> tmparr(v,m);
     318          MaskedArray<Float> tmparr2(tmparr(mask));
     319          Float stdvsq = pow(stddev(tmparr2),2);
     320          stdevsqsum+=1.0/stdvsq;
     321          tmparr /= stdvsq;
     322          outmarr += tmparr;
     323        }
     324        outmarr /= stdevsqsum;
     325        Array<Float> tarr(outmarr.getArray());
     326        Array<Bool> tbarr(outmarr.getMask());
     327        IPosition start(4,i,j,k,0);
     328        IPosition end(4,i,j,k,n-1);
     329        Array<Float> subArr(arr(start,end));
     330        Array<Bool> subMask(barr(start,end));
     331        ArrayAccessor<Float, Axis<0> > aa0(tarr);
     332        ArrayAccessor<Bool, Axis<0> > ba0(tbarr);
     333        ArrayAccessor<Bool, Axis<3> > ba(subMask);
     334        for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
     335          (*aa) = (*aa0);
     336          (*ba) = (*ba0);
     337          aa0++;
     338          ba0++;
     339          ba++;
     340        }
    337341      }
    338342    }
    … …  
    349353}
    350354
    351 CountedPtr<SDMemTable> 
    352 SDMath::averagePol(const CountedPtr<SDMemTable>& in, 
    353                    const Vector<Bool>& mask) {
    354   MaskedArray<Float> marr(in->rowAsMaskedArray(0)); 
     355CountedPtr<SDMemTable>
     356SDMath::averagePol(const CountedPtr<SDMemTable>& in,
     357                   const Vector<Bool>& mask) {
     358  MaskedArray<Float> marr(in->rowAsMaskedArray(0));
    355359  uInt n = in->nChan();
    356360  IPosition ip = marr.shape();
    … …  
    364368      MaskedArray<Float> outmarr(initvec,initmask);
    365369      for (uInt k=0; k<in->nPol();++k) {
    366         IPosition start(4,i,j,k,0);
    367         IPosition end(4,i,j,k,in->nChan()-1);
    368         Array<Float> subArr(arr(start,end));
    369         Array<Bool> subMask(barr(start,end));
    370         Vector<Float> outv;
    371         Vector<Bool> outm;
    372         Vector<Float> v(subArr.nonDegenerate());
    373         Vector<Bool> m(subMask.nonDegenerate());
    374         MaskedArray<Float> tmparr(v,m);
    375         MaskedArray<Float> tmparr2(tmparr(mask));
    376         Float stdvsq = pow(stddev(tmparr2),2);
    377         stdevsqsum+=1.0/stdvsq;
    378         tmparr /= stdvsq;
    379         outmarr += tmparr;
     370        IPosition start(4,i,j,k,0);
     371        IPosition end(4,i,j,k,in->nChan()-1);
     372        Array<Float> subArr(arr(start,end));
     373        Array<Bool> subMask(barr(start,end));
     374        Vector<Float> outv;
     375        Vector<Bool> outm;
     376        Vector<Float> v(subArr.nonDegenerate());
     377        Vector<Bool> m(subMask.nonDegenerate());
     378        MaskedArray<Float> tmparr(v,m);
     379        MaskedArray<Float> tmparr2(tmparr(mask));
     380        Float stdvsq = pow(stddev(tmparr2),2);
     381        stdevsqsum+=1.0/stdvsq;
     382        tmparr /= stdvsq;
     383        outmarr += tmparr;
    380384      }
    381385      outmarr /= stdevsqsum;
    … …  
    384388      // write averaged pol into all pols - fix up to refrom array
    385389      for (uInt k=0; k<in->nPol();++k) {
    386         IPosition start(4,i,j,k,0);
    387         IPosition end(4,i,j,k,n-1);
    388         Array<Float> subArr(arr(start,end));
    389         Array<Bool> subMask(barr(start,end));
    390         ArrayAccessor<Float, Axis<0> > aa0(tarr);       
    391         ArrayAccessor<Bool, Axis<0> > ba0(tbarr);       
    392         ArrayAccessor<Bool, Axis<3> > ba(subMask);     
    393         for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
    394           (*aa) = (*aa0);
    395           (*ba) = (*ba0);
    396           aa0++;
    397           ba0++;
    398           ba++;
    399         }
     390        IPosition start(4,i,j,k,0);
     391        IPosition end(4,i,j,k,n-1);
     392        Array<Float> subArr(arr(start,end));
     393        Array<Bool> subMask(barr(start,end));
     394        ArrayAccessor<Float, Axis<0> > aa0(tarr);
     395        ArrayAccessor<Bool, Axis<0> > ba0(tbarr);
     396        ArrayAccessor<Bool, Axis<3> > ba(subMask);
     397        for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
     398          (*aa) = (*aa0);
     399          (*ba) = (*ba0);
     400          aa0++;
     401          ba0++;
     402          ba++;
     403        }
    400404      }
    401405    }
    … …  
    414418
    415419Float SDMath::rms(const CountedPtr<SDMemTable>& in,
    416                         const std::vector<bool>& mask) {
     420                        const std::vector<bool>& mask) {
    417421  Float rmsval;
    418422  Vector<Bool> msk(mask);
    … …  
    442446}
    443447
    444 CountedPtr<SDMemTable> SDMath::bin(const CountedPtr<SDMemTable>& in,
    445                                           Int width) {
    446  
    447   MaskedArray<Float> marr(in->rowAsMaskedArray(0));
    448   SpectralCoordinate coord(in->getCoordinate(0));
    449   SDContainer sc = in->getSDContainer();
    450   Array<Float> arr = marr.getArray();
    451   Array<Bool> barr = marr.getMask();
    452   SpectralCoordinate outcoord,outcoord2;
    453   MaskedArray<Float> marrout;
    454   ImageUtilities::bin(marrout, outcoord, marr, coord, 3, width);
    455   IPosition ip = marrout.shape();
    456   sc.resize(ip);
    457   sc.putSpectrum(marrout.getArray());
    458   Array<uChar> outflags(ip);
    459   convertArray(outflags,!(marrout.getMask())); 
    460   sc.putFlags(outflags);
    461   MaskedArray<Float> tsys,tsysout;
    462   Array<Bool> dummy(ip);dummy = True;
    463   tsys.setData(sc.getTsys(),dummy);
    464   ImageUtilities::bin(tsysout, outcoord2, marr, outcoord, 3, width);
    465   sc.putTsys(tsysout.getArray());
     448CountedPtr<SDMemTable> SDMath::bin(const CountedPtr<SDMemTable>& in,
     449                                   Int width) {
     450
    466451  SDHeader sh = in->getSDHeader();
     452  SDMemTable* sdmt = new SDMemTable(*in,True);
     453
     454  MaskedArray<Float> tmpmarr(in->rowAsMaskedArray(0));
     455  MaskedArray<Float> tmpmarr2;
     456  SpectralCoordinate outcoord;
     457  IPosition ip;
     458  for (uInt j=0; j<in->nCoordinates(); ++j) {
     459    SpectralCoordinate coord(in->getCoordinate(j));
     460    ImageUtilities::bin(tmpmarr2, outcoord, tmpmarr, coord, 3, width);
     461    ip = tmpmarr2.shape();
     462    sdmt->setCoordinate(outcoord,j);
     463  }
    467464  sh.nchan = ip(3);
    468   SDMemTable* sdmt = new SDMemTable(*in,True);
    469   sdmt->setCoordinate(outcoord,0);
    470   sdmt->putSDContainer(sc);
    471465  sdmt->putSDHeader(sh);
    472   return CountedPtr<SDMemTable>(sdmt);
    473 }
     466
     467  SpectralCoordinate tmpcoord2,tmpcoord;
     468  for (uInt i=0; i < in->nRow(); ++i) {
     469    MaskedArray<Float> marr(in->rowAsMaskedArray(i));
     470    SDContainer sc = in->getSDContainer(i);
     471    Array<Float> arr = marr.getArray();
     472    Array<Bool> barr = marr.getMask();
     473    MaskedArray<Float> marrout;
     474    ImageUtilities::bin(marrout, tmpcoord2, marr, tmpcoord, 3, width);
     475    IPosition ip2 = marrout.shape();
     476    sc.resize(ip2);
     477    sc.putSpectrum(marrout.getArray());
     478    Array<uChar> outflags(ip2);
     479    convertArray(outflags,!(marrout.getMask()));
     480    sc.putFlags(outflags);
     481    MaskedArray<Float> tsys,tsysout;
     482    Array<Bool> dummy(ip);dummy = True;
     483    tsys.setData(sc.getTsys(),dummy);
     484    ImageUtilities::bin(tsysout, tmpcoord2, marr, tmpcoord, 3, width);
     485    sc.putTsys(tsysout.getArray());
     486    sdmt->putSDContainer(sc);
     487  }
     488  return CountedPtr<SDMemTable>(sdmt);
     489}
Note: See TracChangeset for help on using the changeset viewer.