source: trunk/src/SDMath.cc @ 125

Last change on this file since 125 was 125, checked in by mar637, 19 years ago

Moved to casa namespace.
Adjusted the copyright to be ATNF.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.5 KB
Line 
1//#---------------------------------------------------------------------------
2//# SDMath.cc: A collection of single dish mathematical operations
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2004
5//# 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#include <casa/Arrays/VectorIter.h>
45
46#include <tables/Tables/Table.h>
47#include <tables/Tables/ScalarColumn.h>
48#include <tables/Tables/ArrayColumn.h>
49
50#include <scimath/Fitting.h>
51#include <scimath/Fitting/LinearFit.h>
52#include <scimath/Functionals/CompiledFunction.h>
53#include <images/Images/ImageUtilities.h>
54#include <coordinates/Coordinates/SpectralCoordinate.h>
55#include <scimath/Mathematics/AutoDiff.h>
56#include <scimath/Mathematics/AutoDiffMath.h>
57
58#include "MathUtils.h"
59#include "SDContainer.h"
60#include "SDMemTable.h"
61
62#include "SDMath.h"
63
64using namespace casa;
65using namespace asap;
66//using namespace asap::SDMath;
67
68CountedPtr<SDMemTable> SDMath::average(const CountedPtr<SDMemTable>& in) {
69  Table t = in->table();
70  ROArrayColumn<Float> tsys(t, "TSYS");
71  ROScalarColumn<Double> mjd(t, "TIME");
72  ROScalarColumn<String> srcn(t, "SRCNAME");
73  ROScalarColumn<Double> integr(t, "INTERVAL");
74  ROArrayColumn<uInt> freqidc(t, "FREQID");
75  IPosition ip = in->rowAsMaskedArray(0).shape();
76  Array<Float> outarr(ip); outarr =0.0;
77  Array<Float> narr(ip);narr = 0.0;
78  Array<Float> narrinc(ip);narrinc = 1.0;
79
80  Array<Float> tsarr(tsys.shape(0));
81  Array<Float> outtsarr(tsys.shape(0));
82  outtsarr =0.0;
83  tsys.get(0, tsarr);// this is probably unneccessary as tsys should
84  Double tme = 0.0;
85  Double inttime = 0.0;
86
87  for (uInt i=0; i < t.nrow(); i++) {
88    // data stuff
89    MaskedArray<Float> marr(in->rowAsMaskedArray(i));
90    outarr += marr;
91    MaskedArray<Float> n(narrinc,marr.getMask());
92    narr += n;
93    // get
94    tsys.get(i, tsarr);// this is probably unneccessary as tsys should
95    outtsarr += tsarr; // be constant
96    Double tmp;
97    mjd.get(i,tmp);
98    tme += tmp;// average time
99    integr.get(i,tmp);
100    inttime += tmp;
101  }
102  // averaging using mask
103  MaskedArray<Float> nma(narr,(narr > Float(0)));
104  outarr /= nma;
105  Array<Bool> outflagsb = !(nma.getMask());
106  Array<uChar> outflags(outflagsb.shape());
107  convertArray(outflags,outflagsb);
108  SDContainer sc = in->getSDContainer();
109  Int n = t.nrow();
110  outtsarr /= Float(n);
111  sc.timestamp = tme/Double(n);
112  sc.interval = inttime;
113  String tstr; srcn.getScalar(0,tstr);// get sourcename of "mid" point
114  sc.sourcename = tstr;
115  Vector<uInt> tvec;
116  freqidc.get(0,tvec);
117  sc.putFreqMap(tvec);
118  sc.putTsys(outtsarr);
119  sc.scanid = 0;
120  sc.putSpectrum(outarr);
121  sc.putFlags(outflags);
122  SDMemTable* sdmt = new SDMemTable(*in,True);
123  sdmt->putSDContainer(sc);
124  return CountedPtr<SDMemTable>(sdmt);
125}
126
127CountedPtr<SDMemTable>
128SDMath::quotient(const CountedPtr<SDMemTable>& on,
129                 const CountedPtr<SDMemTable>& off) {
130
131  Table ton = on->table();
132  Table toff = off->table();
133  ROArrayColumn<Float> tsys(toff, "TSYS");
134  ROScalarColumn<Double> mjd(ton, "TIME");
135  ROScalarColumn<Double> integr(ton, "INTERVAL");
136  ROScalarColumn<String> srcn(ton, "SRCNAME");
137  ROArrayColumn<uInt> freqidc(ton, "FREQID");
138
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);
160  SDMemTable* sdmt = new SDMemTable(*on, True);
161  sdmt->putSDContainer(sc);
162  return CountedPtr<SDMemTable>(sdmt);
163}
164
165CountedPtr<SDMemTable>
166SDMath::multiply(const CountedPtr<SDMemTable>& in, Float factor) {
167  SDMemTable* sdmt = new SDMemTable(*in);
168  Table t = sdmt->table();
169  ArrayColumn<Float> spec(t,"SPECTRA");
170
171  for (uInt i=0; i < t.nrow(); i++) {
172    // data stuff
173    MaskedArray<Float> marr(sdmt->rowAsMaskedArray(i));
174    marr *= factor;
175    spec.put(i, marr.getArray());
176  }
177  return CountedPtr<SDMemTable>(sdmt);
178}
179
180CountedPtr<SDMemTable>
181SDMath::add(const CountedPtr<SDMemTable>& in, Float offset) {
182  SDMemTable* sdmt = new SDMemTable(*in);
183  Table t = sdmt->table();
184  ArrayColumn<Float> spec(t,"SPECTRA");
185
186  for (uInt i=0; i < t.nrow(); i++) {
187    // data stuff
188    MaskedArray<Float> marr(sdmt->rowAsMaskedArray(i));
189    marr += offset;
190    spec.put(i, marr.getArray());
191  }
192  return CountedPtr<SDMemTable>(sdmt);
193}
194
195
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
215CountedPtr<SDMemTable>
216SDMath::hanning(const CountedPtr<SDMemTable>& in) {
217 
218  //IPosition ip = in->rowAsMaskedArray(0).shape();
219  SDMemTable* sdmt = new SDMemTable(*in,True);
220  for (uInt ri=0; ri < in->nRow(); ++ri) {
221   
222    MaskedArray<Float> marr(in->rowAsMaskedArray(ri));
223
224    Array<Float> arr = marr.getArray();
225    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        }
251    }
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    //
271    Array<uChar> outflags(barr.shape());
272    convertArray(outflags,!barr);
273    SDContainer sc = in->getSDContainer(ri);
274    sc.putSpectrum(arr);
275    sc.putFlags(outflags);
276    sdmt->putSDContainer(sc);
277  }
278  return CountedPtr<SDMemTable>(sdmt);
279}
280
281CountedPtr<SDMemTable>
282SDMath::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;
317        }
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      }
336    }
337  }
338  Array<uChar> outflags(barr.shape());
339  convertArray(outflags,!barr);
340  SDContainer sc = in[0]->getSDContainer();
341  sc.putSpectrum(arr);
342  sc.putFlags(outflags);
343  sc.interval = inttime;
344  SDMemTable* sdmt = new SDMemTable(*in[0],True);
345  sdmt->putSDContainer(sc);
346  return CountedPtr<SDMemTable>(sdmt);
347}
348
349CountedPtr<SDMemTable>
350SDMath::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++;
397        }
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
413Float 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}
441
442CountedPtr<SDMemTable> SDMath::bin(const CountedPtr<SDMemTable>& in,
443                                   Int width) {
444
445  SDHeader sh = in->getSDHeader();
446  SDMemTable* sdmt = new SDMemTable(*in,True);
447
448  MaskedArray<Float> tmpmarr(in->rowAsMaskedArray(0));
449  MaskedArray<Float> tmpmarr2;
450  SpectralCoordinate outcoord;
451  IPosition ip;
452  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);
459  sdmt->putSDHeader(sh);
460
461  SpectralCoordinate tmpcoord2,tmpcoord;
462  for (uInt i=0; i < in->nRow(); ++i) {
463    MaskedArray<Float> marr(in->rowAsMaskedArray(i));
464    SDContainer sc = in->getSDContainer(i);
465    Array<Float> arr = marr.getArray();
466    Array<Bool> barr = marr.getMask();
467    MaskedArray<Float> marrout;
468    ImageUtilities::bin(marrout, tmpcoord2, marr, tmpcoord, 3, width);
469    IPosition ip2 = marrout.shape();
470    sc.resize(ip2);
471    sc.putSpectrum(marrout.getArray());
472    Array<uChar> outflags(ip2);
473    convertArray(outflags,!(marrout.getMask()));
474    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());
480    sdmt->putSDContainer(sc);
481  }
482  return CountedPtr<SDMemTable>(sdmt);
483}
Note: See TracBrowser for help on using the repository browser.