source: trunk/src/SDMath.cc @ 83

Last change on this file since 83 was 83, checked in by mar637, 20 years ago

chnaged namesapce form "atnf_sd" to "asap"

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.6 KB
Line 
1//#---------------------------------------------------------------------------
2//# SDMath.cc: A collection of single dish mathematical operations
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2004
5//# Malte Marquarding, ATNF
6//#
7//# This program is free software; you can redistribute it and/or modify it
8//# under the terms of the GNU General Public License as published by the Free
9//# Software Foundation; either version 2 of the License, or (at your option)
10//# any later version.
11//#
12//# This program is distributed in the hope that it will be useful, but
13//# WITHOUT ANY WARRANTY; without even the implied warranty of
14//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
15//# Public License for more details.
16//#
17//# You should have received a copy of the GNU General Public License along
18//# with this program; if not, write to the Free Software Foundation, Inc.,
19//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20//#
21//# Correspondence concerning this software should be addressed as follows:
22//#        Internet email: Malte.Marquarding@csiro.au
23//#        Postal address: Malte Marquarding,
24//#                        Australia Telescope National Facility,
25//#                        P.O. Box 76,
26//#                        Epping, NSW, 2121,
27//#                        AUSTRALIA
28//#
29//# $Id:
30//#---------------------------------------------------------------------------
31#include <vector>
32
33#include <casa/aips.h>
34#include <casa/BasicSL/String.h>
35#include <casa/Arrays/IPosition.h>
36#include <casa/Arrays/Array.h>
37#include <casa/Arrays/ArrayAccessor.h>
38#include <casa/Arrays/Slice.h>
39#include <casa/Arrays/ArrayMath.h>
40#include <casa/Arrays/ArrayLogical.h>
41#include <casa/Arrays/MaskedArray.h>
42#include <casa/Arrays/MaskArrMath.h>
43#include <casa/Arrays/MaskArrLogi.h>
44
45#include <tables/Tables/Table.h>
46#include <tables/Tables/ScalarColumn.h>
47#include <tables/Tables/ArrayColumn.h>
48
49#include <scimath/Fitting.h>
50#include <scimath/Fitting/LinearFit.h>
51#include <scimath/Functionals/CompiledFunction.h>
52#include <images/Images/ImageUtilities.h>
53#include <coordinates/Coordinates/SpectralCoordinate.h>
54#include <scimath/Mathematics/AutoDiff.h>
55#include <scimath/Mathematics/AutoDiffMath.h>
56
57#include "MathUtils.h"
58#include "SDContainer.h"
59#include "SDMemTable.h"
60
61#include "SDMath.h"
62
63using namespace asap;
64//using namespace asap::SDMath;
65
66CountedPtr<SDMemTable> SDMath::average(const CountedPtr<SDMemTable>& in) {
67  Table t = in->table();
68  ROArrayColumn<Float> tsys(t, "TSYS"); 
69  ROScalarColumn<Double> mjd(t, "TIME");
70  ROScalarColumn<String> srcn(t, "SRCNAME");
71  ROScalarColumn<Double> integr(t, "INTERVAL");
72  ROArrayColumn<uInt> freqidc(t, "FREQID");
73  IPosition ip = in->rowAsMaskedArray(0).shape();
74  Array<Float> outarr(ip); outarr =0.0;
75  Array<Float> narr(ip);narr = 0.0;
76  Array<Float> narrinc(ip);narrinc = 1.0;
77
78  Array<Float> tsarr(tsys.shape(0));
79  Array<Float> outtsarr(tsys.shape(0));
80  outtsarr =0.0;
81  tsys.get(0, tsarr);// this is probably unneccessary as tsys should
82  Double tme = 0.0;
83  Double inttime = 0.0;
84
85  for (uInt i=0; i < t.nrow(); i++) {
86    // data stuff
87    MaskedArray<Float> marr(in->rowAsMaskedArray(i));
88    outarr += marr;
89    MaskedArray<Float> n(narrinc,marr.getMask());
90    narr += n;
91    // get
92    tsys.get(i, tsarr);// this is probably unneccessary as tsys should
93    outtsarr += tsarr; // be constant
94    Double tmp;
95    mjd.get(i,tmp);
96    tme += tmp;// average time
97    integr.get(i,tmp);
98    inttime += tmp;
99  }
100  // averaging using mask
101  MaskedArray<Float> nma(narr,(narr > Float(0)));
102  outarr /= nma;
103  Array<Bool> outflagsb = !(nma.getMask());
104  Array<uChar> outflags(outflagsb.shape());
105  convertArray(outflags,outflagsb);
106  SDContainer sc = in->getSDContainer();
107  Int n = t.nrow();
108  outtsarr /= Float(n);
109  sc.timestamp = tme/Double(n);
110  sc.interval = inttime;
111  String tstr; srcn.getScalar(0,tstr);// get sourcename of "mid" point
112  sc.sourcename = tstr;
113  Vector<uInt> tvec;
114  freqidc.get(0,tvec);
115  sc.putFreqMap(tvec);
116  sc.putTsys(outtsarr);
117  sc.scanid = 0;
118  sc.putSpectrum(outarr);
119  sc.putFlags(outflags); 
120  SDMemTable* sdmt = new SDMemTable(*in,True);
121  sdmt->putSDContainer(sc);
122  return CountedPtr<SDMemTable>(sdmt);
123}
124
125CountedPtr<SDMemTable>
126SDMath::quotient(const CountedPtr<SDMemTable>& on,
127                 const CountedPtr<SDMemTable>& off) {
128 
129  Table ton = on->table();
130  Table toff = off->table();
131  ROArrayColumn<Float> tsys(toff, "TSYS"); 
132  ROScalarColumn<Double> mjd(ton, "TIME");
133  ROScalarColumn<Double> integr(ton, "INTERVAL");
134  ROScalarColumn<String> srcn(ton, "SRCNAME");
135  ROArrayColumn<uInt> freqidc(ton, "FREQID");
136
137  MaskedArray<Float> mon(on->rowAsMaskedArray(0));
138  MaskedArray<Float> moff(off->rowAsMaskedArray(0));
139  IPosition ipon = mon.shape();
140  IPosition ipoff = moff.shape();
141  Array<Float> tsarr;//(tsys.shape(0));
142  tsys.get(0, tsarr);
143  if (ipon != ipoff && ipon != tsarr.shape())
144    cerr << "on/off not conformant" << endl;
145 
146  MaskedArray<Float> tmp = (mon-moff);
147  Array<Float> out(tmp.getArray());
148  out /= moff;
149  out *= tsarr;
150  Array<Bool> outflagsb = !(mon.getMask() && moff.getMask());
151  Array<uChar> outflags(outflagsb.shape());
152  convertArray(outflags,outflagsb);
153  SDContainer sc = on->getSDContainer();
154  sc.putTsys(tsarr);
155  sc.scanid = 0;
156  sc.putSpectrum(out);
157  sc.putFlags(outflags);
158  SDMemTable* sdmt = new SDMemTable(*on, True);
159  sdmt->putSDContainer(sc);
160  return CountedPtr<SDMemTable>(sdmt);
161}
162
163CountedPtr<SDMemTable>
164SDMath::multiply(const CountedPtr<SDMemTable>& in, Float factor) {
165  SDMemTable* sdmt = new SDMemTable(*in);
166  Table t = sdmt->table();
167  ArrayColumn<Float> spec(t,"SPECTRA");
168
169  for (uInt i=0; i < t.nrow(); i++) {
170    // data stuff
171    MaskedArray<Float> marr(sdmt->rowAsMaskedArray(i));
172    marr *= factor;
173    spec.put(i, marr.getArray());
174  }
175  return CountedPtr<SDMemTable>(sdmt);
176}
177
178bool SDMath::fit(Vector<Float>& thefit, const Vector<Float>& data,
179                const Vector<Bool>& mask,
180                const std::string& fitexpr) {
181
182  LinearFit<Float> fitter;
183  Vector<Float> x(data.nelements());
184  indgen(x);
185  CompiledFunction<AutoDiff<Float> > fn;
186  fn.setFunction(String(fitexpr));
187  fitter.setFunction(fn);
188  Vector<Float> out,out1;
189  out = fitter.fit(x,data,&mask);
190  thefit = data;
191  fitter.residual(thefit, x);
192  cout << "Parameter solution = " << out << endl;
193  return True;
194}
195
196CountedPtr<SDMemTable>
197SDMath::baseline(const CountedPtr<SDMemTable>& in,
198                 const std::string& fitexpr,
199                 const std::vector<bool>& mask) {
200 
201  IPosition ip = in->rowAsMaskedArray(0).shape();
202  SDContainer sc = in->getSDContainer();
203  String sname(in->getSourceName());
204  String stim(in->getTime());
205  cout << "Fitting: " << String(fitexpr) << " to "
206       << sname << " [" << stim << "]" << ":" <<endl;
207  MaskedArray<Float> marr(in->rowAsMaskedArray(0));
208  Vector<Bool> inmask(mask);
209  Array<Float> arr = marr.getArray();
210  Array<Bool> barr = marr.getMask();
211  for (uInt i=0; i<in->nBeam();++i) {
212    for (uInt j=0; j<in->nIF();++j) {
213      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        }
228      }
229    }
230  }
231  Array<uChar> outflags(barr.shape());
232  convertArray(outflags,!barr);
233  sc.putSpectrum(arr);
234  sc.putFlags(outflags);
235  SDMemTable* sdmt = new SDMemTable(*in,True);
236  sdmt->putSDContainer(sc);
237  return CountedPtr<SDMemTable>(sdmt);
238}
239
240
241CountedPtr<SDMemTable>
242SDMath::hanning(const CountedPtr<SDMemTable>& in) {
243
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      }
272    }
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
284CountedPtr<SDMemTable>
285SDMath::averages(const Block<CountedPtr<SDMemTable> >& in,
286                 const Vector<Bool>& mask) {
287  IPosition ip = in[0]->rowAsMaskedArray(0).shape();
288  Array<Float> arr(ip);
289  Array<Bool> barr(ip);
290  Double inttime = 0.0;
291 
292  uInt n = in[0]->nChan();
293  for (uInt i=0; i<in[0]->nBeam();++i) {
294    for (uInt j=0; j<in[0]->nIF();++j) {
295      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        }
337      }
338    }
339  }
340  Array<uChar> outflags(barr.shape());
341  convertArray(outflags,!barr);
342  SDContainer sc = in[0]->getSDContainer();
343  sc.putSpectrum(arr);
344  sc.putFlags(outflags);
345  sc.interval = inttime;
346  SDMemTable* sdmt = new SDMemTable(*in[0],True);
347  sdmt->putSDContainer(sc);
348  return CountedPtr<SDMemTable>(sdmt);
349}
350
351CountedPtr<SDMemTable>
352SDMath::averagePol(const CountedPtr<SDMemTable>& in,
353                   const Vector<Bool>& mask) {
354  MaskedArray<Float> marr(in->rowAsMaskedArray(0)); 
355  uInt n = in->nChan();
356  IPosition ip = marr.shape();
357  Array<Float> arr = marr.getArray();
358  Array<Bool> barr = marr.getMask();
359  for (uInt i=0; i<in->nBeam();++i) {
360    for (uInt j=0; j<in->nIF();++j) {
361      Float stdevsqsum = 0.0;
362      Vector<Float> initvec(n);initvec = 0.0;
363      Vector<Bool> initmask(n);initmask = True;
364      MaskedArray<Float> outmarr(initvec,initmask);
365      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;
380      }
381      outmarr /= stdevsqsum;
382      Array<Float> tarr(outmarr.getArray());
383      Array<Bool> tbarr(outmarr.getMask());
384      // write averaged pol into all pols - fix up to refrom array
385      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        }
400      }
401    }
402  }
403
404  Array<uChar> outflags(barr.shape());
405  convertArray(outflags,!barr);
406  SDContainer sc = in->getSDContainer();
407  sc.putSpectrum(arr);
408  sc.putFlags(outflags);
409  SDMemTable* sdmt = new SDMemTable(*in,True);
410  sdmt->putSDContainer(sc);
411  return CountedPtr<SDMemTable>(sdmt);
412}
413
414
415Float SDMath::rms(const CountedPtr<SDMemTable>& in,
416                         const std::vector<bool>& mask) {
417  Float rmsval;
418  Vector<Bool> msk(mask);
419  IPosition ip = in->rowAsMaskedArray(0).shape();
420  MaskedArray<Float> marr(in->rowAsMaskedArray(0));
421
422  Array<Float> arr = marr.getArray();
423  Array<Bool> barr = marr.getMask();
424  uInt i,j,k;
425  i = in->getBeam();
426  j = in->getIF();
427  k = in->getPol();
428  IPosition start(4,i,j,k,0);
429  IPosition end(4,i,j,k,in->nChan()-1);
430  Array<Float> subArr(arr(start,end));
431  Array<Bool> subMask(barr(start,end));
432  Array<Float> v(subArr.nonDegenerate());
433  Array<Bool> m(subMask.nonDegenerate());
434  MaskedArray<Float> tmp;
435  if (msk.nelements() == m.nelements() ) {
436    tmp.setData(v,m&&msk);
437  } else {
438    tmp.setData(v,m);
439  }
440  rmsval = stddev(tmp);
441  return rmsval;
442}
443
444CountedPtr<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());
466  SDHeader sh = in->getSDHeader();
467  sh.nchan = ip(3);
468  SDMemTable* sdmt = new SDMemTable(*in,True);
469  sdmt->setCoordinate(outcoord,0);
470  sdmt->putSDContainer(sc);
471  sdmt->putSDHeader(sh);
472  return CountedPtr<SDMemTable>(sdmt);
473}
Note: See TracBrowser for help on using the repository browser.