source: trunk/src/SDMath.cc @ 107

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

added 'add' function
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.4 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
178CountedPtr<SDMemTable>
179SDMath::add(const CountedPtr<SDMemTable>& in, Float offset) {
180  SDMemTable* sdmt = new SDMemTable(*in);
181  Table t = sdmt->table();
182  ArrayColumn<Float> spec(t,"SPECTRA");
183
184  for (uInt i=0; i < t.nrow(); i++) {
185    // data stuff
186    MaskedArray<Float> marr(sdmt->rowAsMaskedArray(i));
187    marr += offset;
188    spec.put(i, marr.getArray());
189  }
190  return CountedPtr<SDMemTable>(sdmt);
191}
192
193
194
195bool SDMath::fit(Vector<Float>& thefit, const Vector<Float>& data,
196                const Vector<Bool>& mask,
197                const std::string& fitexpr) {
198
199  LinearFit<Float> fitter;
200  Vector<Float> x(data.nelements());
201  indgen(x);
202  CompiledFunction<AutoDiff<Float> > fn;
203  fn.setFunction(String(fitexpr));
204  fitter.setFunction(fn);
205  Vector<Float> out,out1;
206  out = fitter.fit(x,data,&mask);
207  thefit = data;
208  fitter.residual(thefit, x);
209  cout << "Parameter solution = " << out << endl;
210  return True;
211}
212
213CountedPtr<SDMemTable>
214SDMath::baseline(const CountedPtr<SDMemTable>& in,
215                 const std::string& fitexpr,
216                 const std::vector<bool>& mask) {
217
218  IPosition ip = in->rowAsMaskedArray(0).shape();
219  SDContainer sc = in->getSDContainer();
220  String sname(in->getSourceName());
221  String stim(in->getTime());
222  cout << "Fitting: " << String(fitexpr) << " to "
223       << sname << " [" << stim << "]" << ":" <<endl;
224  MaskedArray<Float> marr(in->rowAsMaskedArray(0));
225  Vector<Bool> inmask(mask);
226  Array<Float> arr = marr.getArray();
227  Array<Bool> barr = marr.getMask();
228  for (uInt i=0; i<in->nBeam();++i) {
229    for (uInt j=0; j<in->nIF();++j) {
230      for (uInt k=0; k<in->nPol();++k) {
231        IPosition start(4,i,j,k,0);
232        IPosition end(4,i,j,k,in->nChan()-1);
233        Array<Float> subArr(arr(start,end));
234        Array<Bool> subMask(barr(start,end));
235        Vector<Float> outv;
236        Vector<Float> v(subArr.nonDegenerate());
237        Vector<Bool> m(subMask.nonDegenerate());
238        cout << "\t Polarisation " << k << "\t";
239        SDMath::fit(outv, v, m&&inmask, fitexpr);
240        ArrayAccessor<Float, Axis<0> > aa0(outv);
241        for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
242          (*aa) = (*aa0);
243          aa0++;
244        }
245      }
246    }
247  }
248  Array<uChar> outflags(barr.shape());
249  convertArray(outflags,!barr);
250  sc.putSpectrum(arr);
251  sc.putFlags(outflags);
252  SDMemTable* sdmt = new SDMemTable(*in,True);
253  sdmt->putSDContainer(sc);
254  return CountedPtr<SDMemTable>(sdmt);
255}
256
257
258CountedPtr<SDMemTable>
259SDMath::hanning(const CountedPtr<SDMemTable>& in) {
260
261  //IPosition ip = in->rowAsMaskedArray(0).shape();
262  SDMemTable* sdmt = new SDMemTable(*in,True);
263  for (uInt ri=0; ri < in->nRow(); ++ri) {
264
265    MaskedArray<Float> marr(in->rowAsMaskedArray(ri));
266
267    Array<Float> arr = marr.getArray();
268    Array<Bool> barr = marr.getMask();
269    for (uInt i=0; i<in->nBeam();++i) {
270        for (uInt j=0; j<in->nIF();++j) {
271        for (uInt k=0; k<in->nPol();++k) {
272            IPosition start(4,i,j,k,0);
273            IPosition end(4,i,j,k,in->nChan()-1);
274            Array<Float> subArr(arr(start,end));
275            Array<Bool> subMask(barr(start,end));
276            Vector<Float> outv;
277            Vector<Bool> outm;
278            Vector<Float> v(subArr.nonDegenerate());
279            Vector<Bool> m(subMask.nonDegenerate());
280            ::hanning(outv,outm,v,m);
281            ArrayAccessor<Float, Axis<0> > aa0(outv);
282            ArrayAccessor<Bool, Axis<0> > ba0(outm);
283            ArrayAccessor<Bool, Axis<3> > ba(subMask);
284            for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
285            (*aa) = (*aa0);
286            (*ba) = (*ba0);
287            aa0++;
288            ba0++;
289            ba++;
290            }
291        }
292        }
293    }
294    Array<uChar> outflags(barr.shape());
295    convertArray(outflags,!barr);
296    SDContainer sc = in->getSDContainer(ri);
297    sc.putSpectrum(arr);
298    sc.putFlags(outflags);
299    sdmt->putSDContainer(sc);
300  }
301  return CountedPtr<SDMemTable>(sdmt);
302}
303
304CountedPtr<SDMemTable>
305SDMath::averages(const Block<CountedPtr<SDMemTable> >& in,
306                 const Vector<Bool>& mask) {
307  IPosition ip = in[0]->rowAsMaskedArray(0).shape();
308  Array<Float> arr(ip);
309  Array<Bool> barr(ip);
310  Double inttime = 0.0;
311
312  uInt n = in[0]->nChan();
313  for (uInt i=0; i<in[0]->nBeam();++i) {
314    for (uInt j=0; j<in[0]->nIF();++j) {
315      for (uInt k=0; k<in[0]->nPol();++k) {
316        Float stdevsqsum = 0.0;
317        Vector<Float> initvec(n);initvec = 0.0;
318        Vector<Bool> initmask(n);initmask = True;
319        MaskedArray<Float> outmarr(initvec,initmask);
320        inttime = 0.0;
321        for (uInt bi=0; bi< in.nelements(); ++bi) {
322          MaskedArray<Float> marr(in[bi]->rowAsMaskedArray(0));
323          inttime += in[bi]->getInterval();
324          Array<Float> arr = marr.getArray();
325          Array<Bool> barr = marr.getMask();
326          IPosition start(4,i,j,k,0);
327          IPosition end(4,i,j,k,n-1);
328          Array<Float> subArr(arr(start,end));
329          Array<Bool> subMask(barr(start,end));
330          Vector<Float> outv;
331          Vector<Bool> outm;
332          Vector<Float> v(subArr.nonDegenerate());
333          Vector<Bool> m(subMask.nonDegenerate());
334          MaskedArray<Float> tmparr(v,m);
335          MaskedArray<Float> tmparr2(tmparr(mask));
336          Float stdvsq = pow(stddev(tmparr2),2);
337          stdevsqsum+=1.0/stdvsq;
338          tmparr /= stdvsq;
339          outmarr += tmparr;
340        }
341        outmarr /= stdevsqsum;
342        Array<Float> tarr(outmarr.getArray());
343        Array<Bool> tbarr(outmarr.getMask());
344        IPosition start(4,i,j,k,0);
345        IPosition end(4,i,j,k,n-1);
346        Array<Float> subArr(arr(start,end));
347        Array<Bool> subMask(barr(start,end));
348        ArrayAccessor<Float, Axis<0> > aa0(tarr);
349        ArrayAccessor<Bool, Axis<0> > ba0(tbarr);
350        ArrayAccessor<Bool, Axis<3> > ba(subMask);
351        for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
352          (*aa) = (*aa0);
353          (*ba) = (*ba0);
354          aa0++;
355          ba0++;
356          ba++;
357        }
358      }
359    }
360  }
361  Array<uChar> outflags(barr.shape());
362  convertArray(outflags,!barr);
363  SDContainer sc = in[0]->getSDContainer();
364  sc.putSpectrum(arr);
365  sc.putFlags(outflags);
366  sc.interval = inttime;
367  SDMemTable* sdmt = new SDMemTable(*in[0],True);
368  sdmt->putSDContainer(sc);
369  return CountedPtr<SDMemTable>(sdmt);
370}
371
372CountedPtr<SDMemTable>
373SDMath::averagePol(const CountedPtr<SDMemTable>& in,
374                   const Vector<Bool>& mask) {
375  MaskedArray<Float> marr(in->rowAsMaskedArray(0));
376  uInt n = in->nChan();
377  IPosition ip = marr.shape();
378  Array<Float> arr = marr.getArray();
379  Array<Bool> barr = marr.getMask();
380  for (uInt i=0; i<in->nBeam();++i) {
381    for (uInt j=0; j<in->nIF();++j) {
382      Float stdevsqsum = 0.0;
383      Vector<Float> initvec(n);initvec = 0.0;
384      Vector<Bool> initmask(n);initmask = True;
385      MaskedArray<Float> outmarr(initvec,initmask);
386      for (uInt k=0; k<in->nPol();++k) {
387        IPosition start(4,i,j,k,0);
388        IPosition end(4,i,j,k,in->nChan()-1);
389        Array<Float> subArr(arr(start,end));
390        Array<Bool> subMask(barr(start,end));
391        Vector<Float> outv;
392        Vector<Bool> outm;
393        Vector<Float> v(subArr.nonDegenerate());
394        Vector<Bool> m(subMask.nonDegenerate());
395        MaskedArray<Float> tmparr(v,m);
396        MaskedArray<Float> tmparr2(tmparr(mask));
397        Float stdvsq = pow(stddev(tmparr2),2);
398        stdevsqsum+=1.0/stdvsq;
399        tmparr /= stdvsq;
400        outmarr += tmparr;
401      }
402      outmarr /= stdevsqsum;
403      Array<Float> tarr(outmarr.getArray());
404      Array<Bool> tbarr(outmarr.getMask());
405      // write averaged pol into all pols - fix up to refrom array
406      for (uInt k=0; k<in->nPol();++k) {
407        IPosition start(4,i,j,k,0);
408        IPosition end(4,i,j,k,n-1);
409        Array<Float> subArr(arr(start,end));
410        Array<Bool> subMask(barr(start,end));
411        ArrayAccessor<Float, Axis<0> > aa0(tarr);
412        ArrayAccessor<Bool, Axis<0> > ba0(tbarr);
413        ArrayAccessor<Bool, Axis<3> > ba(subMask);
414        for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
415          (*aa) = (*aa0);
416          (*ba) = (*ba0);
417          aa0++;
418          ba0++;
419          ba++;
420        }
421      }
422    }
423  }
424
425  Array<uChar> outflags(barr.shape());
426  convertArray(outflags,!barr);
427  SDContainer sc = in->getSDContainer();
428  sc.putSpectrum(arr);
429  sc.putFlags(outflags);
430  SDMemTable* sdmt = new SDMemTable(*in,True);
431  sdmt->putSDContainer(sc);
432  return CountedPtr<SDMemTable>(sdmt);
433}
434
435
436Float SDMath::rms(const CountedPtr<SDMemTable>& in,
437                         const std::vector<bool>& mask) {
438  Float rmsval;
439  Vector<Bool> msk(mask);
440  IPosition ip = in->rowAsMaskedArray(0).shape();
441  MaskedArray<Float> marr(in->rowAsMaskedArray(0));
442
443  Array<Float> arr = marr.getArray();
444  Array<Bool> barr = marr.getMask();
445  uInt i,j,k;
446  i = in->getBeam();
447  j = in->getIF();
448  k = in->getPol();
449  IPosition start(4,i,j,k,0);
450  IPosition end(4,i,j,k,in->nChan()-1);
451  Array<Float> subArr(arr(start,end));
452  Array<Bool> subMask(barr(start,end));
453  Array<Float> v(subArr.nonDegenerate());
454  Array<Bool> m(subMask.nonDegenerate());
455  MaskedArray<Float> tmp;
456  if (msk.nelements() == m.nelements() ) {
457    tmp.setData(v,m&&msk);
458  } else {
459    tmp.setData(v,m);
460  }
461  rmsval = stddev(tmp);
462  return rmsval;
463}
464
465CountedPtr<SDMemTable> SDMath::bin(const CountedPtr<SDMemTable>& in,
466                                   Int width) {
467
468  SDHeader sh = in->getSDHeader();
469  SDMemTable* sdmt = new SDMemTable(*in,True);
470
471  MaskedArray<Float> tmpmarr(in->rowAsMaskedArray(0));
472  MaskedArray<Float> tmpmarr2;
473  SpectralCoordinate outcoord;
474  IPosition ip;
475  for (uInt j=0; j<in->nCoordinates(); ++j) {
476    SpectralCoordinate coord(in->getCoordinate(j));
477    ImageUtilities::bin(tmpmarr2, outcoord, tmpmarr, coord, 3, width);
478    ip = tmpmarr2.shape();
479    sdmt->setCoordinate(outcoord,j);
480  }
481  sh.nchan = ip(3);
482  sdmt->putSDHeader(sh);
483
484  SpectralCoordinate tmpcoord2,tmpcoord;
485  for (uInt i=0; i < in->nRow(); ++i) {
486    MaskedArray<Float> marr(in->rowAsMaskedArray(i));
487    SDContainer sc = in->getSDContainer(i);
488    Array<Float> arr = marr.getArray();
489    Array<Bool> barr = marr.getMask();
490    MaskedArray<Float> marrout;
491    ImageUtilities::bin(marrout, tmpcoord2, marr, tmpcoord, 3, width);
492    IPosition ip2 = marrout.shape();
493    sc.resize(ip2);
494    sc.putSpectrum(marrout.getArray());
495    Array<uChar> outflags(ip2);
496    convertArray(outflags,!(marrout.getMask()));
497    sc.putFlags(outflags);
498    MaskedArray<Float> tsys,tsysout;
499    Array<Bool> dummy(ip);dummy = True;
500    tsys.setData(sc.getTsys(),dummy);
501    ImageUtilities::bin(tsysout, tmpcoord2, marr, tmpcoord, 3, width);
502    sc.putTsys(tsysout.getArray());
503    sdmt->putSDContainer(sc);
504  }
505  return CountedPtr<SDMemTable>(sdmt);
506}
Note: See TracBrowser for help on using the repository browser.