source: trunk/src/SDMath.cc@ 95

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

bin function now bins all rows in the table. ALso fixed a couple of bugs.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.0 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
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 SDMemTable* sdmt = new SDMemTable(*in,True);
246 for (uInt ri=0; ri < in->nRow(); ++ri) {
247
248 MaskedArray<Float> marr(in->rowAsMaskedArray(ri));
249
250 Array<Float> arr = marr.getArray();
251 Array<Bool> barr = marr.getMask();
252 for (uInt i=0; i<in->nBeam();++i) {
253 for (uInt j=0; j<in->nIF();++j) {
254 for (uInt k=0; k<in->nPol();++k) {
255 IPosition start(4,i,j,k,0);
256 IPosition end(4,i,j,k,in->nChan()-1);
257 Array<Float> subArr(arr(start,end));
258 Array<Bool> subMask(barr(start,end));
259 Vector<Float> outv;
260 Vector<Bool> outm;
261 Vector<Float> v(subArr.nonDegenerate());
262 Vector<Bool> m(subMask.nonDegenerate());
263 ::hanning(outv,outm,v,m);
264 ArrayAccessor<Float, Axis<0> > aa0(outv);
265 ArrayAccessor<Bool, Axis<0> > ba0(outm);
266 ArrayAccessor<Bool, Axis<3> > ba(subMask);
267 for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
268 (*aa) = (*aa0);
269 (*ba) = (*ba0);
270 aa0++;
271 ba0++;
272 ba++;
273 }
274 }
275 }
276 }
277 Array<uChar> outflags(barr.shape());
278 convertArray(outflags,!barr);
279 SDContainer sc = in->getSDContainer(ri);
280 sc.putSpectrum(arr);
281 sc.putFlags(outflags);
282 sdmt->putSDContainer(sc);
283 }
284 return CountedPtr<SDMemTable>(sdmt);
285}
286
287CountedPtr<SDMemTable>
288SDMath::averages(const Block<CountedPtr<SDMemTable> >& in,
289 const Vector<Bool>& mask) {
290 IPosition ip = in[0]->rowAsMaskedArray(0).shape();
291 Array<Float> arr(ip);
292 Array<Bool> barr(ip);
293 Double inttime = 0.0;
294
295 uInt n = in[0]->nChan();
296 for (uInt i=0; i<in[0]->nBeam();++i) {
297 for (uInt j=0; j<in[0]->nIF();++j) {
298 for (uInt k=0; k<in[0]->nPol();++k) {
299 Float stdevsqsum = 0.0;
300 Vector<Float> initvec(n);initvec = 0.0;
301 Vector<Bool> initmask(n);initmask = True;
302 MaskedArray<Float> outmarr(initvec,initmask);
303 inttime = 0.0;
304 for (uInt bi=0; bi< in.nelements(); ++bi) {
305 MaskedArray<Float> marr(in[bi]->rowAsMaskedArray(0));
306 inttime += in[bi]->getInterval();
307 Array<Float> arr = marr.getArray();
308 Array<Bool> barr = marr.getMask();
309 IPosition start(4,i,j,k,0);
310 IPosition end(4,i,j,k,n-1);
311 Array<Float> subArr(arr(start,end));
312 Array<Bool> subMask(barr(start,end));
313 Vector<Float> outv;
314 Vector<Bool> outm;
315 Vector<Float> v(subArr.nonDegenerate());
316 Vector<Bool> m(subMask.nonDegenerate());
317 MaskedArray<Float> tmparr(v,m);
318 MaskedArray<Float> tmparr2(tmparr(mask));
319 Float stdvsq = pow(stddev(tmparr2),2);
320 stdevsqsum+=1.0/stdvsq;
321 tmparr /= stdvsq;
322 outmarr += tmparr;
323 }
324 outmarr /= stdevsqsum;
325 Array<Float> tarr(outmarr.getArray());
326 Array<Bool> tbarr(outmarr.getMask());
327 IPosition start(4,i,j,k,0);
328 IPosition end(4,i,j,k,n-1);
329 Array<Float> subArr(arr(start,end));
330 Array<Bool> subMask(barr(start,end));
331 ArrayAccessor<Float, Axis<0> > aa0(tarr);
332 ArrayAccessor<Bool, Axis<0> > ba0(tbarr);
333 ArrayAccessor<Bool, Axis<3> > ba(subMask);
334 for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
335 (*aa) = (*aa0);
336 (*ba) = (*ba0);
337 aa0++;
338 ba0++;
339 ba++;
340 }
341 }
342 }
343 }
344 Array<uChar> outflags(barr.shape());
345 convertArray(outflags,!barr);
346 SDContainer sc = in[0]->getSDContainer();
347 sc.putSpectrum(arr);
348 sc.putFlags(outflags);
349 sc.interval = inttime;
350 SDMemTable* sdmt = new SDMemTable(*in[0],True);
351 sdmt->putSDContainer(sc);
352 return CountedPtr<SDMemTable>(sdmt);
353}
354
355CountedPtr<SDMemTable>
356SDMath::averagePol(const CountedPtr<SDMemTable>& in,
357 const Vector<Bool>& mask) {
358 MaskedArray<Float> marr(in->rowAsMaskedArray(0));
359 uInt n = in->nChan();
360 IPosition ip = marr.shape();
361 Array<Float> arr = marr.getArray();
362 Array<Bool> barr = marr.getMask();
363 for (uInt i=0; i<in->nBeam();++i) {
364 for (uInt j=0; j<in->nIF();++j) {
365 Float stdevsqsum = 0.0;
366 Vector<Float> initvec(n);initvec = 0.0;
367 Vector<Bool> initmask(n);initmask = True;
368 MaskedArray<Float> outmarr(initvec,initmask);
369 for (uInt k=0; k<in->nPol();++k) {
370 IPosition start(4,i,j,k,0);
371 IPosition end(4,i,j,k,in->nChan()-1);
372 Array<Float> subArr(arr(start,end));
373 Array<Bool> subMask(barr(start,end));
374 Vector<Float> outv;
375 Vector<Bool> outm;
376 Vector<Float> v(subArr.nonDegenerate());
377 Vector<Bool> m(subMask.nonDegenerate());
378 MaskedArray<Float> tmparr(v,m);
379 MaskedArray<Float> tmparr2(tmparr(mask));
380 Float stdvsq = pow(stddev(tmparr2),2);
381 stdevsqsum+=1.0/stdvsq;
382 tmparr /= stdvsq;
383 outmarr += tmparr;
384 }
385 outmarr /= stdevsqsum;
386 Array<Float> tarr(outmarr.getArray());
387 Array<Bool> tbarr(outmarr.getMask());
388 // write averaged pol into all pols - fix up to refrom array
389 for (uInt k=0; k<in->nPol();++k) {
390 IPosition start(4,i,j,k,0);
391 IPosition end(4,i,j,k,n-1);
392 Array<Float> subArr(arr(start,end));
393 Array<Bool> subMask(barr(start,end));
394 ArrayAccessor<Float, Axis<0> > aa0(tarr);
395 ArrayAccessor<Bool, Axis<0> > ba0(tbarr);
396 ArrayAccessor<Bool, Axis<3> > ba(subMask);
397 for (ArrayAccessor<Float, Axis<3> > aa(subArr); aa != aa.end();++aa) {
398 (*aa) = (*aa0);
399 (*ba) = (*ba0);
400 aa0++;
401 ba0++;
402 ba++;
403 }
404 }
405 }
406 }
407
408 Array<uChar> outflags(barr.shape());
409 convertArray(outflags,!barr);
410 SDContainer sc = in->getSDContainer();
411 sc.putSpectrum(arr);
412 sc.putFlags(outflags);
413 SDMemTable* sdmt = new SDMemTable(*in,True);
414 sdmt->putSDContainer(sc);
415 return CountedPtr<SDMemTable>(sdmt);
416}
417
418
419Float SDMath::rms(const CountedPtr<SDMemTable>& in,
420 const std::vector<bool>& mask) {
421 Float rmsval;
422 Vector<Bool> msk(mask);
423 IPosition ip = in->rowAsMaskedArray(0).shape();
424 MaskedArray<Float> marr(in->rowAsMaskedArray(0));
425
426 Array<Float> arr = marr.getArray();
427 Array<Bool> barr = marr.getMask();
428 uInt i,j,k;
429 i = in->getBeam();
430 j = in->getIF();
431 k = in->getPol();
432 IPosition start(4,i,j,k,0);
433 IPosition end(4,i,j,k,in->nChan()-1);
434 Array<Float> subArr(arr(start,end));
435 Array<Bool> subMask(barr(start,end));
436 Array<Float> v(subArr.nonDegenerate());
437 Array<Bool> m(subMask.nonDegenerate());
438 MaskedArray<Float> tmp;
439 if (msk.nelements() == m.nelements() ) {
440 tmp.setData(v,m&&msk);
441 } else {
442 tmp.setData(v,m);
443 }
444 rmsval = stddev(tmp);
445 return rmsval;
446}
447
448CountedPtr<SDMemTable> SDMath::bin(const CountedPtr<SDMemTable>& in,
449 Int width) {
450
451 SDHeader sh = in->getSDHeader();
452 SDMemTable* sdmt = new SDMemTable(*in,True);
453
454 MaskedArray<Float> tmpmarr(in->rowAsMaskedArray(0));
455 MaskedArray<Float> tmpmarr2;
456 SpectralCoordinate outcoord;
457 IPosition ip;
458 for (uInt j=0; j<in->nCoordinates(); ++j) {
459 SpectralCoordinate coord(in->getCoordinate(j));
460 ImageUtilities::bin(tmpmarr2, outcoord, tmpmarr, coord, 3, width);
461 ip = tmpmarr2.shape();
462 sdmt->setCoordinate(outcoord,j);
463 }
464 sh.nchan = ip(3);
465 sdmt->putSDHeader(sh);
466
467 SpectralCoordinate tmpcoord2,tmpcoord;
468 for (uInt i=0; i < in->nRow(); ++i) {
469 MaskedArray<Float> marr(in->rowAsMaskedArray(i));
470 SDContainer sc = in->getSDContainer(i);
471 Array<Float> arr = marr.getArray();
472 Array<Bool> barr = marr.getMask();
473 MaskedArray<Float> marrout;
474 ImageUtilities::bin(marrout, tmpcoord2, marr, tmpcoord, 3, width);
475 IPosition ip2 = marrout.shape();
476 sc.resize(ip2);
477 sc.putSpectrum(marrout.getArray());
478 Array<uChar> outflags(ip2);
479 convertArray(outflags,!(marrout.getMask()));
480 sc.putFlags(outflags);
481 MaskedArray<Float> tsys,tsysout;
482 Array<Bool> dummy(ip);dummy = True;
483 tsys.setData(sc.getTsys(),dummy);
484 ImageUtilities::bin(tsysout, tmpcoord2, marr, tmpcoord, 3, width);
485 sc.putTsys(tsysout.getArray());
486 sdmt->putSDContainer(sc);
487 }
488 return CountedPtr<SDMemTable>(sdmt);
489}
Note: See TracBrowser for help on using the repository browser.