source: trunk/src/SDMath.cc@ 53

Last change on this file since 53 was 48, checked in by mmarquar, 20 years ago

Fixed various defects. Added averaging of multiple scans, rms, and reworked baseline fitting.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.7 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 <aips/aips.h>
34#include <aips/Utilities/String.h>
35#include <aips/Arrays/IPosition.h>
36#include <aips/Arrays/Array.h>
37#include <aips/Arrays/ArrayAccessor.h>
38#include <aips/Arrays/Slice.h>
39#include <aips/Arrays/ArrayMath.h>
40#include <aips/Arrays/ArrayLogical.h>
41#include <aips/Arrays/MaskedArray.h>
42#include <aips/Arrays/MaskArrMath.h>
43#include <aips/Arrays/MaskArrLogi.h>
44
45#include <aips/Tables/Table.h>
46#include <aips/Tables/ScalarColumn.h>
47#include <aips/Tables/ArrayColumn.h>
48
49#include <aips/Fitting.h>
50#include <trial/Fitting/LinearFit.h>
51#include <trial/Functionals/CompiledFunction.h>
52#include <trial/Images/ImageUtilities.h>
53#include <trial/Coordinates/SpectralCoordinate.h>
54#include <aips/Mathematics/AutoDiff.h>
55#include <aips/Mathematics/AutoDiffMath.h>
56
57#include "MathUtils.h"
58#include "SDContainer.h"
59#include "SDMemTable.h"
60
61#include "SDMath.h"
62
63using namespace atnf_sd;
64
65static CountedPtr<SDMemTable> SDMath::average(const CountedPtr<SDMemTable>& in) {
66 Table t = in->table();
67 ROArrayColumn<Float> tsys(t, "TSYS");
68 ROScalarColumn<Double> mjd(t, "TIME");
69 ROScalarColumn<String> srcn(t, "SRCNAME");
70 ROScalarColumn<Double> integr(t, "INTERVAL");
71 ROArrayColumn<uInt> freqidc(t, "FREQID");
72 IPosition ip = in->rowAsMaskedArray(0).shape();
73 Array<Float> outarr(ip); outarr =0.0;
74 Array<Float> narr(ip);narr = 0.0;
75 Array<Float> narrinc(ip);narrinc = 1.0;
76
77 Array<Float> tsarr(tsys.shape(0));
78 Array<Float> outtsarr(tsys.shape(0));
79 outtsarr =0.0;
80 tsys.get(0, tsarr);// this is probably unneccessary as tsys should
81 Double tme = 0.0;
82 Double inttime = 0.0;
83
84 for (uInt i=0; i < t.nrow(); i++) {
85 // data stuff
86 MaskedArray<Float> marr(in->rowAsMaskedArray(i));
87 outarr += marr;
88 MaskedArray<Float> n(narrinc,marr.getMask());
89 narr += n;
90 // get
91 tsys.get(i, tsarr);// this is probably unneccessary as tsys should
92 outtsarr += tsarr; // be constant
93 Double tmp;
94 mjd.get(i,tmp);
95 tme += tmp;// average time
96 integr.get(i,tmp);
97 inttime += tmp;
98 }
99 // averaging using mask
100 MaskedArray<Float> nma(narr,(narr > Float(0)));
101 outarr /= nma;
102 Array<Bool> outflagsb = !(nma.getMask());
103 Array<uChar> outflags(outflagsb.shape());
104 convertArray(outflags,outflagsb);
105 SDContainer sc = in->getSDContainer();
106 Int n = t.nrow();
107 outtsarr /= Float(n);
108 sc.timestamp = tme/Double(n);
109 sc.interval = inttime;
110 String tstr; srcn.getScalar(0,tstr);// get sourcename of "mid" point
111 sc.sourcename = tstr;
112 Vector<uInt> tvec;
113 freqidc.get(0,tvec);
114 sc.putFreqMap(tvec);
115 sc.putTsys(outtsarr);
116 sc.scanid = 0;
117 sc.putSpectrum(outarr);
118 sc.putFlags(outflags);
119 SDMemTable* sdmt = new SDMemTable(*in,True);
120 sdmt->putSDContainer(sc);
121 return CountedPtr<SDMemTable>(sdmt);
122}
123
124static CountedPtr<SDMemTable>
125SDMath::quotient(const CountedPtr<SDMemTable>& on,
126 const CountedPtr<SDMemTable>& off) {
127
128 Table ton = on->table();
129 Table toff = off->table();
130 ROArrayColumn<Float> tsys(toff, "TSYS");
131 ROScalarColumn<Double> mjd(ton, "TIME");
132 ROScalarColumn<Double> integr(ton, "INTERVAL");
133 ROScalarColumn<String> srcn(ton, "SRCNAME");
134 ROArrayColumn<uInt> freqidc(ton, "FREQID");
135
136 MaskedArray<Float> mon(on->rowAsMaskedArray(0));
137 MaskedArray<Float> moff(off->rowAsMaskedArray(0));
138 IPosition ipon = mon.shape();
139 IPosition ipoff = moff.shape();
140 Array<Float> tsarr;//(tsys.shape(0));
141 tsys.get(0, tsarr);
142 if (ipon != ipoff && ipon != tsarr.shape())
143 cerr << "on/off not conformant" << endl;
144
145 MaskedArray<Float> tmp = (mon-moff);
146 Array<Float> out(tmp.getArray());
147 out /= moff;
148 out *= tsarr;
149 Array<Bool> outflagsb = !(mon.getMask() && moff.getMask());
150 Array<uChar> outflags(outflagsb.shape());
151 convertArray(outflags,outflagsb);
152 SDContainer sc = on->getSDContainer();
153 sc.putTsys(tsarr);
154 sc.scanid = 0;
155 sc.putSpectrum(out);
156 sc.putFlags(outflags);
157 SDMemTable* sdmt = new SDMemTable(*on, True);
158 sdmt->putSDContainer(sc);
159 return CountedPtr<SDMemTable>(sdmt);
160}
161
162static CountedPtr<SDMemTable>
163SDMath::multiply(const CountedPtr<SDMemTable>& in, Float factor) {
164 SDMemTable* sdmt = new SDMemTable(*in);
165 Table t = sdmt->table();
166 ArrayColumn<Float> spec(t,"SPECTRA");
167
168 for (uInt i=0; i < t.nrow(); i++) {
169 // data stuff
170 MaskedArray<Float> marr(sdmt->rowAsMaskedArray(i));
171 marr *= factor;
172 spec.put(i, marr.getArray());
173 }
174 return CountedPtr<SDMemTable>(sdmt);
175}
176
177static bool SDMath::fit(Vector<Float>& thefit, const Vector<Float>& data,
178 const Vector<Bool>& mask,
179 const std::string& fitexpr) {
180
181 LinearFit<Float> fitter;
182 Vector<Float> x(data.nelements());
183 indgen(x);
184 CompiledFunction<AutoDiff<Float> > fn;
185 fn.setFunction(String(fitexpr));
186 fitter.setFunction(fn);
187 Vector<Float> out,out1;
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);
237}
238
239
240static CountedPtr<SDMemTable>
241SDMath::hanning(const CountedPtr<SDMemTable>& in) {
242
243 IPosition ip = in->rowAsMaskedArray(0).shape();
244 MaskedArray<Float> marr(in->rowAsMaskedArray(0));
245
246 Array<Float> arr = marr.getArray();
247 Array<Bool> barr = marr.getMask();
248 for (uInt i=0; i<in->nBeam();++i) {
249 for (uInt j=0; j<in->nIF();++j) {
250 for (uInt k=0; k<in->nPol();++k) {
251 IPosition start(4,i,j,k,0);
252 IPosition end(4,i,j,k,in->nChan()-1);
253 Array<Float> subArr(arr(start,end));
254 Array<Bool> subMask(barr(start,end));
255 Vector<Float> outv;
256 Vector<Bool> outm;
257 Vector<Float> v(subArr.nonDegenerate());
258 Vector<Bool> m(subMask.nonDegenerate());
259 ::hanning(outv,outm,v,m);
260 ArrayAccessor<Float, Axis<0> > aa0(outv);
261 ArrayAccessor<Bool, Axis<0> > ba0(outm);
262 ArrayAccessor<Bool, Axis<3> > ba(subMask);
263 for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
264 (*aa) = (*aa0);
265 (*ba) = (*ba0);
266 aa0++;
267 ba0++;
268 ba++;
269 }
270 }
271 }
272 }
273 Array<uChar> outflags(barr.shape());
274 convertArray(outflags,!barr);
275 SDContainer sc = in->getSDContainer();
276 sc.putSpectrum(arr);
277 sc.putFlags(outflags);
278 SDMemTable* sdmt = new SDMemTable(*in,True);
279 sdmt->putSDContainer(sc);
280 return CountedPtr<SDMemTable>(sdmt);
281}
282
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;
290
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}
Note: See TracBrowser for help on using the repository browser.