Changeset 48 for trunk


Ignore:
Timestamp:
07/15/04 15:48:41 (20 years ago)
Author:
mmarquar
Message:

Fixed various defects. Added averaging of multiple scans, rms, and reworked baseline fitting.

Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/SDMath.cc

    r38 r48  
    5050#include <trial/Fitting/LinearFit.h>
    5151#include <trial/Functionals/CompiledFunction.h>
     52#include <trial/Images/ImageUtilities.h>
     53#include <trial/Coordinates/SpectralCoordinate.h>
    5254#include <aips/Mathematics/AutoDiff.h>
    5355#include <aips/Mathematics/AutoDiffMath.h>
    … …  
    7577  Array<Float> tsarr(tsys.shape(0));
    7678  Array<Float> outtsarr(tsys.shape(0));
     79  outtsarr =0.0;
     80  tsys.get(0, tsarr);// this is probably unneccessary as tsys should
    7781  Double tme = 0.0;
    7882  Double inttime = 0.0;
    … …  
    99103  Array<uChar> outflags(outflagsb.shape());
    100104  convertArray(outflags,outflagsb);
    101   SDContainer sc(ip(0),ip(1),ip(2),ip(3));
    102 
     105  SDContainer sc = in->getSDContainer();
    103106  Int n = t.nrow();
    104   outtsarr /= Float(n/2);
    105   sc.timestamp = tme/Double(n/2);
    106   sc.interval =inttime;
     107  outtsarr /= Float(n);
     108  sc.timestamp = tme/Double(n);
     109  sc.interval = inttime;
    107110  String tstr; srcn.getScalar(0,tstr);// get sourcename of "mid" point
    108111  sc.sourcename = tstr;
    … …  
    110113  freqidc.get(0,tvec);
    111114  sc.putFreqMap(tvec);
     115  sc.putTsys(outtsarr);
    112116  sc.scanid = 0;
    113117  sc.putSpectrum(outarr);
    … …  
    134138  IPosition ipon = mon.shape();
    135139  IPosition ipoff = moff.shape();
    136   Array<Float> tsarr(tsys.shape(0));
    137 
     140  Array<Float> tsarr;//(tsys.shape(0));
     141  tsys.get(0, tsarr);
    138142  if (ipon != ipoff && ipon != tsarr.shape())
    139143    cerr << "on/off not conformant" << endl;
    140144 
    141   //IPosition test = mon.shape()/2;
    142145  MaskedArray<Float> tmp = (mon-moff);
    143146  Array<Float> out(tmp.getArray());
    … …  
    147150  Array<uChar> outflags(outflagsb.shape());
    148151  convertArray(outflags,outflagsb);
    149 
    150   SDContainer sc(ipon(0),ipon(1),ipon(2),ipon(3));
    151   String tstr; srcn.getScalar(0,tstr);// get sourcename of "on" scan
    152   sc.sourcename = tstr;
    153   Double tme; mjd.getScalar(0,tme);// get time of "on" scan
    154   sc.timestamp = tme;
    155   integr.getScalar(0,tme);
    156   sc.interval = tme;
    157   Vector<uInt> tvec;
    158   freqidc.get(0,tvec);
    159   sc.putFreqMap(tvec);
     152  SDContainer sc = on->getSDContainer();
     153  sc.putTsys(tsarr);
    160154  sc.scanid = 0;
    161155  sc.putSpectrum(out);
    162   sc.putFlags(outflags); 
     156  sc.putFlags(outflags);
    163157  SDMemTable* sdmt = new SDMemTable(*on, True);
    164158  sdmt->putSDContainer(sc);
    165159  return CountedPtr<SDMemTable>(sdmt);
    166  
    167 }
     160}
     161
    168162static CountedPtr<SDMemTable>
    169163SDMath::multiply(const CountedPtr<SDMemTable>& in, Float factor) {
    … …  
    180174  return CountedPtr<SDMemTable>(sdmt);
    181175}
    182 static std::vector<float> SDMath::baseline(const CountedPtr<SDMemTable>& in,
    183                                            const std::string& fitexpr) {
    184   cout << "Fitting: " << fitexpr << endl;
    185   Table t = in->table();
     176
     177static bool SDMath::fit(Vector<Float>& thefit, const Vector<Float>& data,
     178                      const Vector<Bool>& mask,
     179                      const std::string& fitexpr) {
     180
    186181  LinearFit<Float> fitter;
    187   Vector<Float> y;
    188   in->getSpectrum(y, 0);
    189   Vector<Bool> m;
    190   in->getMask(m, 0);
    191   Vector<Float> x(y.nelements());
     182  Vector<Float> x(data.nelements());
    192183  indgen(x);
    193184  CompiledFunction<AutoDiff<Float> > fn;
    … …  
    195186  fitter.setFunction(fn);
    196187  Vector<Float> out,out1;
    197   out = fitter.fit(x,y,&m);
    198   out1 = y;
    199   fitter.residual(out1,x);
    200   cout << "solution =" << out << endl;
    201   std::vector<float> fitted;
    202   out1.tovector(fitted);
    203   return fitted;
     188  out = fitter.fit(x,data,&mask);
     189  thefit = data;
     190  fitter.residual(thefit, x);
     191  cout << "Parameter solution = " << out << endl;
     192  return True;
     193}
     194
     195static CountedPtr<SDMemTable>
     196SDMath::baseline(const CountedPtr<SDMemTable>& in,
     197                 const std::string& fitexpr,
     198                 const std::vector<bool>& mask) {
     199 
     200  IPosition ip = in->rowAsMaskedArray(0).shape();
     201  SDContainer sc = in->getSDContainer();
     202  String sname(in->getSourceName());
     203  String stim(in->getTime());
     204  cout << "Fitting: " << String(fitexpr) << " to "
     205       << sname << " [" << stim << "]" << ":" <<endl;
     206  MaskedArray<Float> marr(in->rowAsMaskedArray(0));
     207  Vector<Bool> inmask(mask);
     208  Array<Float> arr = marr.getArray();
     209  Array<Bool> barr = marr.getMask();
     210  for (uInt i=0; i<in->nBeam();++i) {
     211    for (uInt j=0; j<in->nIF();++j) {
     212      for (uInt k=0; k<in->nPol();++k) {
     213        IPosition start(4,i,j,k,0);
     214        IPosition end(4,i,j,k,in->nChan()-1);
     215        Array<Float> subArr(arr(start,end));
     216        Array<Bool> subMask(barr(start,end));
     217        Vector<Float> outv;
     218        Vector<Float> v(subArr.nonDegenerate());
     219        Vector<Bool> m(subMask.nonDegenerate());
     220        cout << "\t Polarisation " << k << "\t";
     221        SDMath::fit(outv, v, m&&inmask, fitexpr);
     222        ArrayAccessor<Float, Axis<0> > aa0(outv);
     223        for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
     224          (*aa) = (*aa0);
     225          aa0++;
     226        }
     227      }
     228    }
     229  }
     230  Array<uChar> outflags(barr.shape());
     231  convertArray(outflags,!barr);
     232  sc.putSpectrum(arr);
     233  sc.putFlags(outflags);
     234  SDMemTable* sdmt = new SDMemTable(*in,True);
     235  sdmt->putSDContainer(sc);
     236  return CountedPtr<SDMemTable>(sdmt);
    204237}
    205238
    … …  
    207240static CountedPtr<SDMemTable>
    208241SDMath::hanning(const CountedPtr<SDMemTable>& in) {
     242
    209243  IPosition ip = in->rowAsMaskedArray(0).shape();
    210244  MaskedArray<Float> marr(in->rowAsMaskedArray(0));
    … …  
    247281}
    248282
    249 /*
    250 static Float SDMath::rms(const SDMemTable& in, uInt whichRow) {
    251   Table t = in.table(); 
    252   MaskedArray<Float> marr(in.rowAsMaskedArray(whichRow,True));
     283static CountedPtr<SDMemTable>
     284SDMath::averages(const Block<CountedPtr<SDMemTable> >& in,
     285                 const Vector<Bool>& mask) {
     286  IPosition ip = in[0]->rowAsMaskedArray(0).shape();
     287  Array<Float> arr(ip);
     288  Array<Bool> barr(ip);
     289  Double inttime = 0.0;
    253290 
    254 }
    255 */
     291  uInt n = in[0]->nChan();
     292  for (uInt i=0; i<in[0]->nBeam();++i) {
     293    for (uInt j=0; j<in[0]->nIF();++j) {
     294      for (uInt k=0; k<in[0]->nPol();++k) {
     295        Float stdevsqsum = 0.0;
     296        Vector<Float> initvec(n);initvec = 0.0;
     297        Vector<Bool> initmask(n);initmask = True;
     298        MaskedArray<Float> outmarr(initvec,initmask);
     299        for (uInt bi=0; bi< in.nelements(); ++bi) {
     300          MaskedArray<Float> marr(in[bi]->rowAsMaskedArray(0));
     301          inttime += in[bi]->getInterval();
     302          Array<Float> arr = marr.getArray();
     303          Array<Bool> barr = marr.getMask();
     304          IPosition start(4,i,j,k,0);
     305          IPosition end(4,i,j,k,n-1);
     306          Array<Float> subArr(arr(start,end));
     307          Array<Bool> subMask(barr(start,end));
     308          Vector<Float> outv;
     309          Vector<Bool> outm;
     310          Vector<Float> v(subArr.nonDegenerate());
     311          Vector<Bool> m(subMask.nonDegenerate());
     312          MaskedArray<Float> tmparr(v,m);
     313          MaskedArray<Float> tmparr2(tmparr(mask));
     314          Float stdvsq = pow(stddev(tmparr2),2);
     315          stdevsqsum+=1.0/stdvsq;
     316          tmparr /= stdvsq;
     317          outmarr += tmparr;
     318        }
     319        outmarr /= stdevsqsum;
     320        Array<Float> tarr(outmarr.getArray());
     321        Array<Bool> tbarr(outmarr.getMask());
     322        IPosition start(4,i,j,k,0);
     323        IPosition end(4,i,j,k,n-1);
     324        Array<Float> subArr(arr(start,end));
     325        Array<Bool> subMask(barr(start,end));
     326        ArrayAccessor<Float, Axis<0> > aa0(tarr);       
     327        ArrayAccessor<Bool, Axis<0> > ba0(tbarr);       
     328        ArrayAccessor<Bool, Axis<3> > ba(subMask);     
     329        for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
     330          (*aa) = (*aa0);
     331          (*ba) = (*ba0);
     332          aa0++;
     333          ba0++;
     334          ba++;
     335        }
     336      }
     337    }
     338  }
     339  Array<uChar> outflags(barr.shape());
     340  convertArray(outflags,!barr);
     341  SDContainer sc = in[0]->getSDContainer();
     342  sc.putSpectrum(arr);
     343  sc.putFlags(outflags);
     344  sc.interval = inttime;
     345  SDMemTable* sdmt = new SDMemTable(*in[0],True);
     346  sdmt->putSDContainer(sc);
     347  return CountedPtr<SDMemTable>(sdmt);
     348}
     349
     350static CountedPtr<SDMemTable>
     351SDMath::averagePol(const CountedPtr<SDMemTable>& in,
     352                   const Vector<Bool>& mask) {
     353  MaskedArray<Float> marr(in->rowAsMaskedArray(0)); 
     354  uInt n = in->nChan();
     355  IPosition ip = marr.shape();
     356  Array<Float> arr = marr.getArray();
     357  Array<Bool> barr = marr.getMask();
     358  for (uInt i=0; i<in->nBeam();++i) {
     359    for (uInt j=0; j<in->nIF();++j) {
     360      Float stdevsqsum = 0.0;
     361      Vector<Float> initvec(n);initvec = 0.0;
     362      Vector<Bool> initmask(n);initmask = True;
     363      MaskedArray<Float> outmarr(initvec,initmask);
     364      for (uInt k=0; k<in->nPol();++k) {
     365        IPosition start(4,i,j,k,0);
     366        IPosition end(4,i,j,k,in->nChan()-1);
     367        Array<Float> subArr(arr(start,end));
     368        Array<Bool> subMask(barr(start,end));
     369        Vector<Float> outv;
     370        Vector<Bool> outm;
     371        Vector<Float> v(subArr.nonDegenerate());
     372        Vector<Bool> m(subMask.nonDegenerate());
     373        MaskedArray<Float> tmparr(v,m);
     374        MaskedArray<Float> tmparr2(tmparr(mask));
     375        Float stdvsq = pow(stddev(tmparr2),2);
     376        stdevsqsum+=1.0/stdvsq;
     377        tmparr /= stdvsq;
     378        outmarr += tmparr;
     379      }
     380      outmarr /= stdevsqsum;
     381      Array<Float> tarr(outmarr.getArray());
     382      Array<Bool> tbarr(outmarr.getMask());
     383      // write averaged pol into all pols - fix up to refrom array
     384      for (uInt k=0; k<in->nPol();++k) {
     385        IPosition start(4,i,j,k,0);
     386        IPosition end(4,i,j,k,n-1);
     387        Array<Float> subArr(arr(start,end));
     388        Array<Bool> subMask(barr(start,end));
     389        ArrayAccessor<Float, Axis<0> > aa0(tarr);       
     390        ArrayAccessor<Bool, Axis<0> > ba0(tbarr);       
     391        ArrayAccessor<Bool, Axis<3> > ba(subMask);     
     392        for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
     393          (*aa) = (*aa0);
     394          (*ba) = (*ba0);
     395          aa0++;
     396          ba0++;
     397          ba++;
     398        }
     399      }
     400    }
     401  }
     402
     403  Array<uChar> outflags(barr.shape());
     404  convertArray(outflags,!barr);
     405  SDContainer sc = in->getSDContainer();
     406  sc.putSpectrum(arr);
     407  sc.putFlags(outflags);
     408  SDMemTable* sdmt = new SDMemTable(*in,True);
     409  sdmt->putSDContainer(sc);
     410  return CountedPtr<SDMemTable>(sdmt);
     411}
     412
     413
     414static Float SDMath::rms(const CountedPtr<SDMemTable>& in,
     415                         const std::vector<bool>& mask) {
     416  Float rmsval;
     417  Vector<Bool> msk(mask);
     418  IPosition ip = in->rowAsMaskedArray(0).shape();
     419  MaskedArray<Float> marr(in->rowAsMaskedArray(0));
     420
     421  Array<Float> arr = marr.getArray();
     422  Array<Bool> barr = marr.getMask();
     423  uInt i,j,k;
     424  i = in->getBeam();
     425  j = in->getIF();
     426  k = in->getPol();
     427  IPosition start(4,i,j,k,0);
     428  IPosition end(4,i,j,k,in->nChan()-1);
     429  Array<Float> subArr(arr(start,end));
     430  Array<Bool> subMask(barr(start,end));
     431  Array<Float> v(subArr.nonDegenerate());
     432  Array<Bool> m(subMask.nonDegenerate());
     433  MaskedArray<Float> tmp;
     434  if (msk.nelements() == m.nelements() ) {
     435    tmp.setData(v,m&&msk);
     436  } else {
     437    tmp.setData(v,m);
     438  }
     439  rmsval = stddev(tmp);
     440  return rmsval;
     441}
     442
     443static CountedPtr<SDMemTable> SDMath::bin(const CountedPtr<SDMemTable>& in,
     444                                          Int width) {
     445 
     446  MaskedArray<Float> marr(in->rowAsMaskedArray(0));
     447  SpectralCoordinate coord(in->getCoordinate(0));
     448  SDContainer sc = in->getSDContainer();
     449  Array<Float> arr = marr.getArray();
     450  Array<Bool> barr = marr.getMask();
     451  SpectralCoordinate outcoord,outcoord2;
     452  MaskedArray<Float> marrout;
     453  ImageUtilities::bin(marrout, outcoord, marr, coord, 3, width);
     454  IPosition ip = marrout.shape();
     455  sc.resize(ip);
     456  sc.putSpectrum(marrout.getArray());
     457  Array<uChar> outflags(ip);
     458  convertArray(outflags,!(marrout.getMask())); 
     459  sc.putFlags(outflags);
     460  MaskedArray<Float> tsys,tsysout;
     461  Array<Bool> dummy(ip);dummy = True;
     462  tsys.setData(sc.getTsys(),dummy);
     463  ImageUtilities::bin(tsysout, outcoord2, marr, outcoord, 3, width);
     464  sc.putTsys(tsysout.getArray());
     465  SDHeader sh = in->getSDHeader();
     466  sh.nchan = ip(3);
     467  SDMemTable* sdmt = new SDMemTable(*in,True);
     468  sdmt->setCoordinate(outcoord,0);
     469  sdmt->putSDContainer(sc);
     470  sdmt->putSDHeader(sh);
     471  return CountedPtr<SDMemTable>(sdmt);
     472}

  • trunk/src/SDMath.h

    r38 r48  
    4848                                         Float factor);
    4949 
    50   static std::vector<float> baseline(const CountedPtr<SDMemTable>& in,
    51                                      const std::string& fitexpr );
     50  static CountedPtr<SDMemTable> baseline(const CountedPtr<SDMemTable>& in,
     51                                         const std::string& fitexpr,
     52                                         const std::vector<bool>& mask);
    5253  static CountedPtr<SDMemTable> hanning(const CountedPtr<SDMemTable>& in);
     54
     55  static CountedPtr<SDMemTable>
     56  averages(const Block<CountedPtr<SDMemTable> >& in,
     57           const Vector<Bool>& mask);
     58
     59  static CountedPtr<SDMemTable>
     60  averagePol(const CountedPtr<SDMemTable>& in, const Vector<Bool>& mask);
     61
     62  static Float rms(const CountedPtr<SDMemTable>& in,
     63                   const std::vector<bool>& mask);
     64 
     65  static CountedPtr<SDMemTable> bin(const CountedPtr<SDMemTable>& in,
     66                                    Int width);
     67 
     68private:
     69  static bool fit(Vector<Float>& thefit, const Vector<Float>& data,
     70                  const Vector<Bool>& mask, const std::string& fitexpr);
    5371 
    5472};
Note: See TracChangeset for help on using the changeset viewer.