source: trunk/src/SDMath.cc@ 71

Last change on this file since 71 was 66, checked in by mmarquar, 20 years ago

Replaced class and statics with a namespace. Functions now global in the namespace SDMath
  • 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 <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//using namespace atnf_sd::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.