source: trunk/src/SDMath.cc@ 128

Last change on this file since 128 was 125, checked in by mar637, 21 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.