source: trunk/src/SDMath.cc@ 113

Last change on this file since 113 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.