source: trunk/src/Scantable.cpp@ 2873

Last change on this file since 2873 was 2869, checked in by Takeshi Nakazato, 11 years ago

New Development: No

JIRA Issue: Yes CAS-5841

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: test_sdsave, test_importasdm_sd, all sd regressions

Put in Release Notes: Yes

Module(s): sd

Description: Describe your changes here...

Filler/writer are changed so that SCANNO is consistent with
SCAN_NUMBER in MS and scanNumber in ASDM.

It affects several regressions and unit tests so that tests
are necessary to be updated. This change should be verified
by the updated unit tests/regression tests.


  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 156.8 KB
Line 
1//
2// C++ Implementation: Scantable
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2005-2013
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12#include <map>
13#include <sys/time.h>
14
15#include <atnf/PKSIO/SrcType.h>
16
17#include <casa/aips.h>
18#include <casa/iomanip.h>
19#include <casa/iostream.h>
20#include <casa/OS/File.h>
21#include <casa/OS/Path.h>
22#include <casa/Logging/LogIO.h>
23#include <casa/Arrays/Array.h>
24#include <casa/Arrays/ArrayAccessor.h>
25#include <casa/Arrays/ArrayLogical.h>
26#include <casa/Arrays/ArrayMath.h>
27#include <casa/Arrays/MaskArrMath.h>
28#include <casa/Arrays/Slice.h>
29#include <casa/Arrays/Vector.h>
30#include <casa/Arrays/VectorSTLIterator.h>
31#include <casa/BasicMath/Math.h>
32#include <casa/BasicSL/Constants.h>
33#include <casa/Containers/RecordField.h>
34#include <casa/Logging/LogIO.h>
35#include <casa/Quanta/MVAngle.h>
36#include <casa/Quanta/MVTime.h>
37#include <casa/Utilities/GenSort.h>
38
39#include <coordinates/Coordinates/CoordinateUtil.h>
40
41// needed to avoid error in .tcc
42#include <measures/Measures/MCDirection.h>
43//
44#include <measures/Measures/MDirection.h>
45#include <measures/Measures/MEpoch.h>
46#include <measures/Measures/MFrequency.h>
47#include <measures/Measures/MeasRef.h>
48#include <measures/Measures/MeasTable.h>
49#include <measures/TableMeasures/ScalarMeasColumn.h>
50#include <measures/TableMeasures/TableMeasDesc.h>
51#include <measures/TableMeasures/TableMeasRefDesc.h>
52#include <measures/TableMeasures/TableMeasValueDesc.h>
53
54#include <tables/Tables/ArrColDesc.h>
55#include <tables/Tables/ExprNode.h>
56#include <tables/Tables/ScaColDesc.h>
57#include <tables/Tables/SetupNewTab.h>
58#include <tables/Tables/TableCopy.h>
59#include <tables/Tables/TableDesc.h>
60#include <tables/Tables/TableIter.h>
61#include <tables/Tables/TableParse.h>
62#include <tables/Tables/TableRecord.h>
63#include <tables/Tables/TableRow.h>
64#include <tables/Tables/TableVector.h>
65
66#include "MathUtils.h"
67#include "STAttr.h"
68#include "STBaselineTable.h"
69#include "STLineFinder.h"
70#include "STPolCircular.h"
71#include "STPolLinear.h"
72#include "STPolStokes.h"
73#include "STUpgrade.h"
74#include "STFitter.h"
75#include "Scantable.h"
76
77#define debug 1
78
79using namespace casa;
80
81namespace asap {
82
83std::map<std::string, STPol::STPolFactory *> Scantable::factories_;
84
85void Scantable::initFactories() {
86 if ( factories_.empty() ) {
87 Scantable::factories_["linear"] = &STPolLinear::myFactory;
88 Scantable::factories_["circular"] = &STPolCircular::myFactory;
89 Scantable::factories_["stokes"] = &STPolStokes::myFactory;
90 }
91}
92
93Scantable::Scantable(Table::TableType ttype) :
94 type_(ttype)
95{
96 initFactories();
97 setupMainTable();
98 freqTable_ = STFrequencies(*this);
99 table_.rwKeywordSet().defineTable("FREQUENCIES", freqTable_.table());
100 weatherTable_ = STWeather(*this);
101 table_.rwKeywordSet().defineTable("WEATHER", weatherTable_.table());
102 focusTable_ = STFocus(*this);
103 table_.rwKeywordSet().defineTable("FOCUS", focusTable_.table());
104 tcalTable_ = STTcal(*this);
105 table_.rwKeywordSet().defineTable("TCAL", tcalTable_.table());
106 moleculeTable_ = STMolecules(*this);
107 table_.rwKeywordSet().defineTable("MOLECULES", moleculeTable_.table());
108 historyTable_ = STHistory(*this);
109 table_.rwKeywordSet().defineTable("HISTORY", historyTable_.table());
110 fitTable_ = STFit(*this);
111 table_.rwKeywordSet().defineTable("FIT", fitTable_.table());
112 table_.tableInfo().setType( "Scantable" ) ;
113 originalTable_ = table_;
114 attach();
115}
116
117Scantable::Scantable(const std::string& name, Table::TableType ttype) :
118 type_(ttype)
119{
120 initFactories();
121
122 Table tab(name, Table::Update);
123 uInt version = tab.keywordSet().asuInt("VERSION");
124 if (version != version_) {
125 STUpgrade upgrader(version_);
126 LogIO os( LogOrigin( "Scantable" ) ) ;
127 os << LogIO::WARN
128 << name << " data format version " << version
129 << " is deprecated" << endl
130 << "Running upgrade."<< endl
131 << LogIO::POST ;
132 std::string outname = upgrader.upgrade(name);
133 if ( outname != name ) {
134 os << LogIO::WARN
135 << "Data will be loaded from " << outname << " instead of "
136 << name << LogIO::POST ;
137 tab = Table(outname, Table::Update ) ;
138 }
139 }
140 if ( type_ == Table::Memory ) {
141 table_ = tab.copyToMemoryTable(generateName());
142 } else {
143 table_ = tab;
144 }
145 table_.tableInfo().setType( "Scantable" ) ;
146
147 attachSubtables();
148 originalTable_ = table_;
149 attach();
150}
151/*
152Scantable::Scantable(const std::string& name, Table::TableType ttype) :
153 type_(ttype)
154{
155 initFactories();
156 Table tab(name, Table::Update);
157 uInt version = tab.keywordSet().asuInt("VERSION");
158 if (version != version_) {
159 throw(AipsError("Unsupported version of ASAP file."));
160 }
161 if ( type_ == Table::Memory ) {
162 table_ = tab.copyToMemoryTable(generateName());
163 } else {
164 table_ = tab;
165 }
166
167 attachSubtables();
168 originalTable_ = table_;
169 attach();
170}
171*/
172
173Scantable::Scantable( const Scantable& other, bool clear )
174{
175 // with or without data
176 String newname = String(generateName());
177 type_ = other.table_.tableType();
178 if ( other.table_.tableType() == Table::Memory ) {
179 if ( clear ) {
180 table_ = TableCopy::makeEmptyMemoryTable(newname,
181 other.table_, True);
182 } else {
183 table_ = other.table_.copyToMemoryTable(newname);
184 }
185 } else {
186 other.table_.deepCopy(newname, Table::New, False,
187 other.table_.endianFormat(),
188 Bool(clear));
189 table_ = Table(newname, Table::Update);
190 table_.markForDelete();
191 }
192 table_.tableInfo().setType( "Scantable" ) ;
193 /// @todo reindex SCANNO, recompute nbeam, nif, npol
194 if ( clear ) copySubtables(other);
195 attachSubtables();
196 originalTable_ = table_;
197 attach();
198}
199
200void Scantable::copySubtables(const Scantable& other) {
201 Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
202 TableCopy::copyRows(t, other.freqTable_.table());
203 t = table_.rwKeywordSet().asTable("FOCUS");
204 TableCopy::copyRows(t, other.focusTable_.table());
205 t = table_.rwKeywordSet().asTable("WEATHER");
206 TableCopy::copyRows(t, other.weatherTable_.table());
207 t = table_.rwKeywordSet().asTable("TCAL");
208 TableCopy::copyRows(t, other.tcalTable_.table());
209 t = table_.rwKeywordSet().asTable("MOLECULES");
210 TableCopy::copyRows(t, other.moleculeTable_.table());
211 t = table_.rwKeywordSet().asTable("HISTORY");
212 TableCopy::copyRows(t, other.historyTable_.table());
213 t = table_.rwKeywordSet().asTable("FIT");
214 TableCopy::copyRows(t, other.fitTable_.table());
215}
216
217void Scantable::attachSubtables()
218{
219 freqTable_ = STFrequencies(table_);
220 focusTable_ = STFocus(table_);
221 weatherTable_ = STWeather(table_);
222 tcalTable_ = STTcal(table_);
223 moleculeTable_ = STMolecules(table_);
224 historyTable_ = STHistory(table_);
225 fitTable_ = STFit(table_);
226}
227
228Scantable::~Scantable()
229{
230}
231
232void Scantable::setupMainTable()
233{
234 TableDesc td("", "1", TableDesc::Scratch);
235 td.comment() = "An ASAP Scantable";
236 td.rwKeywordSet().define("VERSION", uInt(version_));
237
238 // n Cycles
239 td.addColumn(ScalarColumnDesc<uInt>("SCANNO"));
240 // new index every nBeam x nIF x nPol
241 td.addColumn(ScalarColumnDesc<uInt>("CYCLENO"));
242
243 td.addColumn(ScalarColumnDesc<uInt>("BEAMNO"));
244 td.addColumn(ScalarColumnDesc<uInt>("IFNO"));
245 // linear, circular, stokes
246 td.rwKeywordSet().define("POLTYPE", String("linear"));
247 td.addColumn(ScalarColumnDesc<uInt>("POLNO"));
248
249 td.addColumn(ScalarColumnDesc<uInt>("FREQ_ID"));
250 td.addColumn(ScalarColumnDesc<uInt>("MOLECULE_ID"));
251
252 ScalarColumnDesc<Int> refbeamnoColumn("REFBEAMNO");
253 refbeamnoColumn.setDefault(Int(-1));
254 td.addColumn(refbeamnoColumn);
255
256 ScalarColumnDesc<uInt> flagrowColumn("FLAGROW");
257 flagrowColumn.setDefault(uInt(0));
258 td.addColumn(flagrowColumn);
259
260 td.addColumn(ScalarColumnDesc<Double>("TIME"));
261 TableMeasRefDesc measRef(MEpoch::UTC); // UTC as default
262 TableMeasValueDesc measVal(td, "TIME");
263 TableMeasDesc<MEpoch> mepochCol(measVal, measRef);
264 mepochCol.write(td);
265
266 td.addColumn(ScalarColumnDesc<Double>("INTERVAL"));
267
268 td.addColumn(ScalarColumnDesc<String>("SRCNAME"));
269 // Type of source (on=0, off=1, other=-1)
270 ScalarColumnDesc<Int> stypeColumn("SRCTYPE");
271 stypeColumn.setDefault(Int(-1));
272 td.addColumn(stypeColumn);
273 td.addColumn(ScalarColumnDesc<String>("FIELDNAME"));
274
275 //The actual Data Vectors
276 td.addColumn(ArrayColumnDesc<Float>("SPECTRA"));
277 td.addColumn(ArrayColumnDesc<uChar>("FLAGTRA"));
278 td.addColumn(ArrayColumnDesc<Float>("TSYS"));
279
280 td.addColumn(ArrayColumnDesc<Double>("DIRECTION",
281 IPosition(1,2),
282 ColumnDesc::Direct));
283 TableMeasRefDesc mdirRef(MDirection::J2000); // default
284 TableMeasValueDesc tmvdMDir(td, "DIRECTION");
285 // the TableMeasDesc gives the column a type
286 TableMeasDesc<MDirection> mdirCol(tmvdMDir, mdirRef);
287 // a uder set table type e.g. GALCTIC, B1950 ...
288 td.rwKeywordSet().define("DIRECTIONREF", String("J2000"));
289 // writing create the measure column
290 mdirCol.write(td);
291 td.addColumn(ScalarColumnDesc<Float>("AZIMUTH"));
292 td.addColumn(ScalarColumnDesc<Float>("ELEVATION"));
293 td.addColumn(ScalarColumnDesc<Float>("OPACITY"));
294
295 td.addColumn(ScalarColumnDesc<uInt>("TCAL_ID"));
296 ScalarColumnDesc<Int> fitColumn("FIT_ID");
297 fitColumn.setDefault(Int(-1));
298 td.addColumn(fitColumn);
299
300 td.addColumn(ScalarColumnDesc<uInt>("FOCUS_ID"));
301 td.addColumn(ScalarColumnDesc<uInt>("WEATHER_ID"));
302
303 // columns which just get dragged along, as they aren't used in asap
304 td.addColumn(ScalarColumnDesc<Double>("SRCVELOCITY"));
305 td.addColumn(ArrayColumnDesc<Double>("SRCPROPERMOTION"));
306 td.addColumn(ArrayColumnDesc<Double>("SRCDIRECTION"));
307 td.addColumn(ArrayColumnDesc<Double>("SCANRATE"));
308
309 td.rwKeywordSet().define("OBSMODE", String(""));
310
311 // Now create Table SetUp from the description.
312 SetupNewTable aNewTab(generateName(), td, Table::Scratch);
313 table_ = Table(aNewTab, type_, 0);
314 originalTable_ = table_;
315}
316
317void Scantable::attach()
318{
319 timeCol_.attach(table_, "TIME");
320 srcnCol_.attach(table_, "SRCNAME");
321 srctCol_.attach(table_, "SRCTYPE");
322 specCol_.attach(table_, "SPECTRA");
323 flagsCol_.attach(table_, "FLAGTRA");
324 tsysCol_.attach(table_, "TSYS");
325 cycleCol_.attach(table_,"CYCLENO");
326 scanCol_.attach(table_, "SCANNO");
327 beamCol_.attach(table_, "BEAMNO");
328 ifCol_.attach(table_, "IFNO");
329 polCol_.attach(table_, "POLNO");
330 integrCol_.attach(table_, "INTERVAL");
331 azCol_.attach(table_, "AZIMUTH");
332 elCol_.attach(table_, "ELEVATION");
333 dirCol_.attach(table_, "DIRECTION");
334 fldnCol_.attach(table_, "FIELDNAME");
335 rbeamCol_.attach(table_, "REFBEAMNO");
336
337 mweatheridCol_.attach(table_,"WEATHER_ID");
338 mfitidCol_.attach(table_,"FIT_ID");
339 mfreqidCol_.attach(table_, "FREQ_ID");
340 mtcalidCol_.attach(table_, "TCAL_ID");
341 mfocusidCol_.attach(table_, "FOCUS_ID");
342 mmolidCol_.attach(table_, "MOLECULE_ID");
343
344 //Add auxiliary column for row-based flagging (CAS-1433 Wataru Kawasaki)
345 attachAuxColumnDef(flagrowCol_, "FLAGROW", 0);
346
347}
348
349template<class T, class T2>
350void Scantable::attachAuxColumnDef(ScalarColumn<T>& col,
351 const String& colName,
352 const T2& defValue)
353{
354 try {
355 col.attach(table_, colName);
356 } catch (TableError& err) {
357 String errMesg = err.getMesg();
358 if (errMesg == "Table column " + colName + " is unknown") {
359 table_.addColumn(ScalarColumnDesc<T>(colName));
360 col.attach(table_, colName);
361 col.fillColumn(static_cast<T>(defValue));
362 } else {
363 throw;
364 }
365 } catch (...) {
366 throw;
367 }
368}
369
370template<class T, class T2>
371void Scantable::attachAuxColumnDef(ArrayColumn<T>& col,
372 const String& colName,
373 const Array<T2>& defValue)
374{
375 try {
376 col.attach(table_, colName);
377 } catch (TableError& err) {
378 String errMesg = err.getMesg();
379 if (errMesg == "Table column " + colName + " is unknown") {
380 table_.addColumn(ArrayColumnDesc<T>(colName));
381 col.attach(table_, colName);
382
383 int size = 0;
384 ArrayIterator<T2>& it = defValue.begin();
385 while (it != defValue.end()) {
386 ++size;
387 ++it;
388 }
389 IPosition ip(1, size);
390 Array<T>& arr(ip);
391 for (int i = 0; i < size; ++i)
392 arr[i] = static_cast<T>(defValue[i]);
393
394 col.fillColumn(arr);
395 } else {
396 throw;
397 }
398 } catch (...) {
399 throw;
400 }
401}
402
403void Scantable::setHeader(const STHeader& sdh)
404{
405 table_.rwKeywordSet().define("nIF", sdh.nif);
406 table_.rwKeywordSet().define("nBeam", sdh.nbeam);
407 table_.rwKeywordSet().define("nPol", sdh.npol);
408 table_.rwKeywordSet().define("nChan", sdh.nchan);
409 table_.rwKeywordSet().define("Observer", sdh.observer);
410 table_.rwKeywordSet().define("Project", sdh.project);
411 table_.rwKeywordSet().define("Obstype", sdh.obstype);
412 table_.rwKeywordSet().define("AntennaName", sdh.antennaname);
413 table_.rwKeywordSet().define("AntennaPosition", sdh.antennaposition);
414 table_.rwKeywordSet().define("Equinox", sdh.equinox);
415 table_.rwKeywordSet().define("FreqRefFrame", sdh.freqref);
416 table_.rwKeywordSet().define("FreqRefVal", sdh.reffreq);
417 table_.rwKeywordSet().define("Bandwidth", sdh.bandwidth);
418 table_.rwKeywordSet().define("UTC", sdh.utc);
419 table_.rwKeywordSet().define("FluxUnit", sdh.fluxunit);
420 table_.rwKeywordSet().define("Epoch", sdh.epoch);
421 table_.rwKeywordSet().define("POLTYPE", sdh.poltype);
422}
423
424STHeader Scantable::getHeader() const
425{
426 STHeader sdh;
427 table_.keywordSet().get("nBeam",sdh.nbeam);
428 table_.keywordSet().get("nIF",sdh.nif);
429 table_.keywordSet().get("nPol",sdh.npol);
430 table_.keywordSet().get("nChan",sdh.nchan);
431 table_.keywordSet().get("Observer", sdh.observer);
432 table_.keywordSet().get("Project", sdh.project);
433 table_.keywordSet().get("Obstype", sdh.obstype);
434 table_.keywordSet().get("AntennaName", sdh.antennaname);
435 table_.keywordSet().get("AntennaPosition", sdh.antennaposition);
436 table_.keywordSet().get("Equinox", sdh.equinox);
437 table_.keywordSet().get("FreqRefFrame", sdh.freqref);
438 table_.keywordSet().get("FreqRefVal", sdh.reffreq);
439 table_.keywordSet().get("Bandwidth", sdh.bandwidth);
440 table_.keywordSet().get("UTC", sdh.utc);
441 table_.keywordSet().get("FluxUnit", sdh.fluxunit);
442 table_.keywordSet().get("Epoch", sdh.epoch);
443 table_.keywordSet().get("POLTYPE", sdh.poltype);
444 return sdh;
445}
446
447void Scantable::setSourceType( int stype )
448{
449 if ( stype < 0 || stype > 1 )
450 throw(AipsError("Illegal sourcetype."));
451 TableVector<Int> tabvec(table_, "SRCTYPE");
452 tabvec = Int(stype);
453}
454
455void Scantable::setSourceName( const std::string& name )
456{
457 TableVector<String> tabvec(table_, "SRCNAME");
458 tabvec = name;
459}
460
461bool Scantable::conformant( const Scantable& other )
462{
463 return this->getHeader().conformant(other.getHeader());
464}
465
466
467
468std::string Scantable::formatSec(Double x) const
469{
470 Double xcop = x;
471 MVTime mvt(xcop/24./3600.); // make days
472
473 if (x < 59.95)
474 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_HM, 7)+"s";
475 else if (x < 3599.95)
476 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_H,7)+" ";
477 else {
478 ostringstream oss;
479 oss << setw(2) << std::right << setprecision(1) << mvt.hour();
480 oss << ":" << mvt.string(MVTime::TIME_CLEAN_NO_H,7) << " ";
481 return String(oss);
482 }
483};
484
485std::string Scantable::formatDirection(const MDirection& md) const
486{
487 Vector<Double> t = md.getAngle(Unit(String("rad"))).getValue();
488 Int prec = 7;
489
490 String ref = md.getRefString();
491 MVAngle mvLon(t[0]);
492 String sLon = mvLon.string(MVAngle::TIME,prec);
493 uInt tp = md.getRef().getType();
494 if (tp == MDirection::GALACTIC ||
495 tp == MDirection::SUPERGAL ) {
496 sLon = mvLon(0.0).string(MVAngle::ANGLE_CLEAN,prec);
497 }
498 MVAngle mvLat(t[1]);
499 String sLat = mvLat.string(MVAngle::ANGLE+MVAngle::DIG2,prec);
500 return ref + String(" ") + sLon + String(" ") + sLat;
501}
502
503
504std::string Scantable::getFluxUnit() const
505{
506 return table_.keywordSet().asString("FluxUnit");
507}
508
509void Scantable::setFluxUnit(const std::string& unit)
510{
511 String tmp(unit);
512 Unit tU(tmp);
513 if (tU==Unit("K") || tU==Unit("Jy")) {
514 table_.rwKeywordSet().define(String("FluxUnit"), tmp);
515 } else {
516 throw AipsError("Illegal unit - must be compatible with Jy or K");
517 }
518}
519
520void Scantable::setInstrument(const std::string& name)
521{
522 bool throwIt = true;
523 // create an Instrument to see if this is valid
524 STAttr::convertInstrument(name, throwIt);
525 String nameU(name);
526 nameU.upcase();
527 table_.rwKeywordSet().define(String("AntennaName"), nameU);
528}
529
530void Scantable::setFeedType(const std::string& feedtype)
531{
532 if ( Scantable::factories_.find(feedtype) == Scantable::factories_.end() ) {
533 std::string msg = "Illegal feed type "+ feedtype;
534 throw(casa::AipsError(msg));
535 }
536 table_.rwKeywordSet().define(String("POLTYPE"), feedtype);
537}
538
539MPosition Scantable::getAntennaPosition() const
540{
541 Vector<Double> antpos;
542 table_.keywordSet().get("AntennaPosition", antpos);
543 MVPosition mvpos(antpos(0),antpos(1),antpos(2));
544 return MPosition(mvpos);
545}
546
547void Scantable::makePersistent(const std::string& filename)
548{
549 String inname(filename);
550 Path path(inname);
551 /// @todo reindex SCANNO, recompute nbeam, nif, npol
552 inname = path.expandedName();
553 // 2011/03/04 TN
554 // We can comment out this workaround since the essential bug is
555 // fixed in casacore (r20889 in google code).
556 table_.deepCopy(inname, Table::New);
557// // WORKAROUND !!! for Table bug
558// // Remove when fixed in casacore
559// if ( table_.tableType() == Table::Memory && !selector_.empty() ) {
560// Table tab = table_.copyToMemoryTable(generateName());
561// tab.deepCopy(inname, Table::New);
562// tab.markForDelete();
563//
564// } else {
565// table_.deepCopy(inname, Table::New);
566// }
567}
568
569int Scantable::nbeam( int scanno ) const
570{
571 if ( scanno < 0 ) {
572 Int n;
573 table_.keywordSet().get("nBeam",n);
574 return int(n);
575 } else {
576 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
577 Table t = table_(table_.col("SCANNO") == scanno);
578 ROTableRow row(t);
579 const TableRecord& rec = row.get(0);
580 Table subt = t( t.col("IFNO") == Int(rec.asuInt("IFNO"))
581 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
582 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
583 ROTableVector<uInt> v(subt, "BEAMNO");
584 return int(v.nelements());
585 }
586 return 0;
587}
588
589int Scantable::nif( int scanno ) const
590{
591 if ( scanno < 0 ) {
592 Int n;
593 table_.keywordSet().get("nIF",n);
594 return int(n);
595 } else {
596 // take the first POLNO,BEAMNO,CYCLENO as nbeam shouldn't vary with these
597 Table t = table_(table_.col("SCANNO") == scanno);
598 ROTableRow row(t);
599 const TableRecord& rec = row.get(0);
600 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
601 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
602 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
603 if ( subt.nrow() == 0 ) return 0;
604 ROTableVector<uInt> v(subt, "IFNO");
605 return int(v.nelements());
606 }
607 return 0;
608}
609
610int Scantable::npol( int scanno ) const
611{
612 if ( scanno < 0 ) {
613 Int n;
614 table_.keywordSet().get("nPol",n);
615 return n;
616 } else {
617 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
618 Table t = table_(table_.col("SCANNO") == scanno);
619 ROTableRow row(t);
620 const TableRecord& rec = row.get(0);
621 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
622 && t.col("IFNO") == Int(rec.asuInt("IFNO"))
623 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
624 if ( subt.nrow() == 0 ) return 0;
625 ROTableVector<uInt> v(subt, "POLNO");
626 return int(v.nelements());
627 }
628 return 0;
629}
630
631int Scantable::ncycle( int scanno ) const
632{
633 if ( scanno < 0 ) {
634 Block<String> cols(2);
635 cols[0] = "SCANNO";
636 cols[1] = "CYCLENO";
637 TableIterator it(table_, cols);
638 int n = 0;
639 while ( !it.pastEnd() ) {
640 ++n;
641 ++it;
642 }
643 return n;
644 } else {
645 Table t = table_(table_.col("SCANNO") == scanno);
646 ROTableRow row(t);
647 const TableRecord& rec = row.get(0);
648 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
649 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
650 && t.col("IFNO") == Int(rec.asuInt("IFNO")) );
651 if ( subt.nrow() == 0 ) return 0;
652 return int(subt.nrow());
653 }
654 return 0;
655}
656
657
658int Scantable::nrow( int scanno ) const
659{
660 return int(table_.nrow());
661}
662
663int Scantable::nchan( int ifno ) const
664{
665 if ( ifno < 0 ) {
666 Int n;
667 table_.keywordSet().get("nChan",n);
668 return int(n);
669 } else {
670 // take the first SCANNO,POLNO,BEAMNO,CYCLENO as nbeam shouldn't
671 // vary with these
672 Table t = table_(table_.col("IFNO") == ifno, 1);
673 if ( t.nrow() == 0 ) return 0;
674 ROArrayColumn<Float> v(t, "SPECTRA");
675 return v.shape(0)(0);
676 }
677 return 0;
678}
679
680int Scantable::nscan() const {
681 Vector<uInt> scannos(scanCol_.getColumn());
682 uInt nout = genSort( scannos, Sort::Ascending,
683 Sort::QuickSort|Sort::NoDuplicates );
684 return int(nout);
685}
686
687int Scantable::getChannels(int whichrow) const
688{
689 return specCol_.shape(whichrow)(0);
690}
691
692int Scantable::getBeam(int whichrow) const
693{
694 return beamCol_(whichrow);
695}
696
697std::vector<uint> Scantable::getNumbers(const ScalarColumn<uInt>& col) const
698{
699 Vector<uInt> nos(col.getColumn());
700 uInt n = genSort( nos, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates );
701 nos.resize(n, True);
702 std::vector<uint> stlout;
703 nos.tovector(stlout);
704 return stlout;
705}
706
707int Scantable::getIF(int whichrow) const
708{
709 return ifCol_(whichrow);
710}
711
712int Scantable::getPol(int whichrow) const
713{
714 return polCol_(whichrow);
715}
716
717std::string Scantable::formatTime(const MEpoch& me, bool showdate) const
718{
719 return formatTime(me, showdate, 0);
720}
721
722std::string Scantable::formatTime(const MEpoch& me, bool showdate, uInt prec) const
723{
724 MVTime mvt(me.getValue());
725 if (showdate)
726 //mvt.setFormat(MVTime::YMD);
727 mvt.setFormat(MVTime::YMD, prec);
728 else
729 //mvt.setFormat(MVTime::TIME);
730 mvt.setFormat(MVTime::TIME, prec);
731 ostringstream oss;
732 oss << mvt;
733 return String(oss);
734}
735
736void Scantable::calculateAZEL()
737{
738 LogIO os( LogOrigin( "Scantable", "calculateAZEL()", WHERE ) ) ;
739 MPosition mp = getAntennaPosition();
740 MEpoch::ROScalarColumn timeCol(table_, "TIME");
741 ostringstream oss;
742 oss << mp;
743 os << "Computed azimuth/elevation using " << endl
744 << String(oss) << endl;
745 for (Int i=0; i<nrow(); ++i) {
746 MEpoch me = timeCol(i);
747 MDirection md = getDirection(i);
748 os << " Time: " << formatTime(me,False)
749 << " Direction: " << formatDirection(md)
750 << endl << " => ";
751 MeasFrame frame(mp, me);
752 Vector<Double> azel =
753 MDirection::Convert(md, MDirection::Ref(MDirection::AZEL,
754 frame)
755 )().getAngle("rad").getValue();
756 azCol_.put(i,Float(azel[0]));
757 elCol_.put(i,Float(azel[1]));
758 os << "azel: " << azel[0]/C::pi*180.0 << " "
759 << azel[1]/C::pi*180.0 << " (deg)" << LogIO::POST;
760 }
761}
762
763void Scantable::clip(const Float uthres, const Float dthres, bool clipoutside, bool unflag)
764{
765 for (uInt i=0; i<table_.nrow(); ++i) {
766 Vector<uChar> flgs = flagsCol_(i);
767 srchChannelsToClip(i, uthres, dthres, clipoutside, unflag, flgs);
768 flagsCol_.put(i, flgs);
769 }
770}
771
772std::vector<bool> Scantable::getClipMask(int whichrow, const Float uthres, const Float dthres, bool clipoutside, bool unflag)
773{
774 Vector<uChar> flags;
775 flagsCol_.get(uInt(whichrow), flags);
776 srchChannelsToClip(uInt(whichrow), uthres, dthres, clipoutside, unflag, flags);
777 Vector<Bool> bflag(flags.shape());
778 convertArray(bflag, flags);
779 //bflag = !bflag;
780
781 std::vector<bool> mask;
782 bflag.tovector(mask);
783 return mask;
784}
785
786void Scantable::srchChannelsToClip(uInt whichrow, const Float uthres, const Float dthres, bool clipoutside, bool unflag,
787 Vector<uChar> flgs)
788{
789 Vector<Float> spcs = specCol_(whichrow);
790 uInt nchannel = spcs.nelements();
791 if (spcs.nelements() != nchannel) {
792 throw(AipsError("Data has incorrect number of channels"));
793 }
794 uChar userflag = 1 << 7;
795 if (unflag) {
796 userflag = 0 << 7;
797 }
798 if (clipoutside) {
799 for (uInt j = 0; j < nchannel; ++j) {
800 Float spc = spcs(j);
801 if ((spc >= uthres) || (spc <= dthres)) {
802 flgs(j) = userflag;
803 }
804 }
805 } else {
806 for (uInt j = 0; j < nchannel; ++j) {
807 Float spc = spcs(j);
808 if ((spc < uthres) && (spc > dthres)) {
809 flgs(j) = userflag;
810 }
811 }
812 }
813}
814
815
816void Scantable::flag( int whichrow, const std::vector<bool>& msk, bool unflag ) {
817 std::vector<bool>::const_iterator it;
818 uInt ntrue = 0;
819 if (whichrow >= int(table_.nrow()) ) {
820 throw(AipsError("Invalid row number"));
821 }
822 for (it = msk.begin(); it != msk.end(); ++it) {
823 if ( *it ) {
824 ntrue++;
825 }
826 }
827 //if ( selector_.empty() && (msk.size() == 0 || msk.size() == ntrue) )
828 if ( whichrow == -1 && !unflag && selector_.empty() && (msk.size() == 0 || msk.size() == ntrue) )
829 throw(AipsError("Trying to flag whole scantable."));
830 uChar userflag = 1 << 7;
831 if ( unflag ) {
832 userflag = 0 << 7;
833 }
834 if (whichrow > -1 ) {
835 applyChanFlag(uInt(whichrow), msk, userflag);
836 } else {
837 for ( uInt i=0; i<table_.nrow(); ++i) {
838 applyChanFlag(i, msk, userflag);
839 }
840 }
841}
842
843void Scantable::applyChanFlag( uInt whichrow, const std::vector<bool>& msk, uChar flagval )
844{
845 if (whichrow >= table_.nrow() ) {
846 throw( casa::indexError<int>( whichrow, "asap::Scantable::applyChanFlag: Invalid row number" ) );
847 }
848 Vector<uChar> flgs = flagsCol_(whichrow);
849 if ( msk.size() == 0 ) {
850 flgs = flagval;
851 flagsCol_.put(whichrow, flgs);
852 return;
853 }
854 if ( int(msk.size()) != nchan( getIF(whichrow) ) ) {
855 throw(AipsError("Mask has incorrect number of channels."));
856 }
857 if ( flgs.nelements() != msk.size() ) {
858 throw(AipsError("Mask has incorrect number of channels."
859 " Probably varying with IF. Please flag per IF"));
860 }
861 std::vector<bool>::const_iterator it;
862 uInt j = 0;
863 for (it = msk.begin(); it != msk.end(); ++it) {
864 if ( *it ) {
865 flgs(j) = flagval;
866 }
867 ++j;
868 }
869 flagsCol_.put(whichrow, flgs);
870}
871
872void Scantable::flagRow(const std::vector<uInt>& rows, bool unflag)
873{
874 if (selector_.empty() && (rows.size() == table_.nrow()) && !unflag)
875 throw(AipsError("Trying to flag whole scantable."));
876
877 uInt rowflag = (unflag ? 0 : 1);
878 std::vector<uInt>::const_iterator it;
879 for (it = rows.begin(); it != rows.end(); ++it)
880 flagrowCol_.put(*it, rowflag);
881}
882
883std::vector<bool> Scantable::getMask(int whichrow) const
884{
885 Vector<uChar> flags;
886 flagsCol_.get(uInt(whichrow), flags);
887 Vector<Bool> bflag(flags.shape());
888 convertArray(bflag, flags);
889 bflag = !bflag;
890 std::vector<bool> mask;
891 bflag.tovector(mask);
892 return mask;
893}
894
895std::vector<float> Scantable::getSpectrum( int whichrow,
896 const std::string& poltype ) const
897{
898 LogIO os( LogOrigin( "Scantable", "getSpectrum()", WHERE ) ) ;
899
900 String ptype = poltype;
901 if (poltype == "" ) ptype = getPolType();
902 if ( whichrow < 0 || whichrow >= nrow() )
903 throw(AipsError("Illegal row number."));
904 std::vector<float> out;
905 Vector<Float> arr;
906 uInt requestedpol = polCol_(whichrow);
907 String basetype = getPolType();
908 if ( ptype == basetype ) {
909 specCol_.get(whichrow, arr);
910 } else {
911 CountedPtr<STPol> stpol(STPol::getPolClass(Scantable::factories_,
912 basetype));
913 uInt row = uInt(whichrow);
914 stpol->setSpectra(getPolMatrix(row));
915 Float fang,fhand;
916 fang = focusTable_.getTotalAngle(mfocusidCol_(row));
917 fhand = focusTable_.getFeedHand(mfocusidCol_(row));
918 stpol->setPhaseCorrections(fang, fhand);
919 arr = stpol->getSpectrum(requestedpol, ptype);
920 }
921 if ( arr.nelements() == 0 )
922
923 os << "Not enough polarisations present to do the conversion."
924 << LogIO::POST;
925 arr.tovector(out);
926 return out;
927}
928
929void Scantable::setSpectrum( const std::vector<float>& spec,
930 int whichrow )
931{
932 Vector<Float> spectrum(spec);
933 Vector<Float> arr;
934 specCol_.get(whichrow, arr);
935 if ( spectrum.nelements() != arr.nelements() )
936 throw AipsError("The spectrum has incorrect number of channels.");
937 specCol_.put(whichrow, spectrum);
938}
939
940
941String Scantable::generateName()
942{
943 return (File::newUniqueName("./","temp")).baseName();
944}
945
946const casa::Table& Scantable::table( ) const
947{
948 return table_;
949}
950
951casa::Table& Scantable::table( )
952{
953 return table_;
954}
955
956std::string Scantable::getPolType() const
957{
958 return table_.keywordSet().asString("POLTYPE");
959}
960
961void Scantable::unsetSelection()
962{
963 table_ = originalTable_;
964 attach();
965 selector_.reset();
966}
967
968void Scantable::setSelection( const STSelector& selection )
969{
970 Table tab = const_cast<STSelector&>(selection).apply(originalTable_);
971 if ( tab.nrow() == 0 ) {
972 throw(AipsError("Selection contains no data. Not applying it."));
973 }
974 table_ = tab;
975 attach();
976// tab.rwKeywordSet().define("nBeam",(Int)(getBeamNos().size())) ;
977// vector<uint> selectedIFs = getIFNos() ;
978// Int newnIF = selectedIFs.size() ;
979// tab.rwKeywordSet().define("nIF",newnIF) ;
980// if ( newnIF != 0 ) {
981// Int newnChan = 0 ;
982// for ( Int i = 0 ; i < newnIF ; i++ ) {
983// Int nChan = nchan( selectedIFs[i] ) ;
984// if ( newnChan > nChan )
985// newnChan = nChan ;
986// }
987// tab.rwKeywordSet().define("nChan",newnChan) ;
988// }
989// tab.rwKeywordSet().define("nPol",(Int)(getPolNos().size())) ;
990 selector_ = selection;
991}
992
993
994std::string Scantable::headerSummary()
995{
996 // Format header info
997// STHeader sdh;
998// sdh = getHeader();
999// sdh.print();
1000 ostringstream oss;
1001 oss.flags(std::ios_base::left);
1002 String tmp;
1003 // Project
1004 table_.keywordSet().get("Project", tmp);
1005 oss << setw(15) << "Project:" << tmp << endl;
1006 // Observation date
1007 oss << setw(15) << "Obs Date:" << getTime(-1,true) << endl;
1008 // Observer
1009 oss << setw(15) << "Observer:"
1010 << table_.keywordSet().asString("Observer") << endl;
1011 // Antenna Name
1012 table_.keywordSet().get("AntennaName", tmp);
1013 oss << setw(15) << "Antenna Name:" << tmp << endl;
1014 // Obs type
1015 table_.keywordSet().get("Obstype", tmp);
1016 // Records (nrow)
1017 oss << setw(15) << "Data Records:" << table_.nrow() << " rows" << endl;
1018 oss << setw(15) << "Obs. Type:" << tmp << endl;
1019 // Beams, IFs, Polarizations, and Channels
1020 oss << setw(15) << "Beams:" << setw(4) << nbeam() << endl
1021 << setw(15) << "IFs:" << setw(4) << nif() << endl
1022 << setw(15) << "Polarisations:" << setw(4) << npol()
1023 << "(" << getPolType() << ")" << endl
1024 << setw(15) << "Channels:" << nchan() << endl;
1025 // Flux unit
1026 table_.keywordSet().get("FluxUnit", tmp);
1027 oss << setw(15) << "Flux Unit:" << tmp << endl;
1028 // Abscissa Unit
1029 oss << setw(15) << "Abscissa:" << getAbcissaLabel(0) << endl;
1030 // Selection
1031 oss << selector_.print() << endl;
1032
1033 return String(oss);
1034}
1035
1036void Scantable::summary( const std::string& filename )
1037{
1038 ostringstream oss;
1039 ofstream ofs;
1040 LogIO ols(LogOrigin("Scantable", "summary", WHERE));
1041
1042 if (filename != "")
1043 ofs.open( filename.c_str(), ios::out );
1044
1045 oss << endl;
1046 oss << asap::SEPERATOR << endl;
1047 oss << " Scan Table Summary" << endl;
1048 oss << asap::SEPERATOR << endl;
1049
1050 // Format header info
1051 oss << headerSummary();
1052 oss << endl;
1053
1054 if (table_.nrow() <= 0){
1055 oss << asap::SEPERATOR << endl;
1056 oss << "The MAIN table is empty: there are no data!!!" << endl;
1057 oss << asap::SEPERATOR << endl;
1058
1059 ols << String(oss) << LogIO::POST;
1060 if (ofs) {
1061 ofs << String(oss) << flush;
1062 ofs.close();
1063 }
1064 return;
1065 }
1066
1067
1068
1069 // main table
1070 String dirtype = "Position ("
1071 + getDirectionRefString()
1072 + ")";
1073 oss.flags(std::ios_base::left);
1074 oss << setw(5) << "Scan"
1075 << setw(15) << "Source"
1076 << setw(35) << "Time range"
1077 << setw(2) << "" << setw(7) << "Int[s]"
1078 << setw(7) << "Record"
1079 << setw(8) << "SrcType"
1080 << setw(8) << "FreqIDs"
1081 << setw(7) << "MolIDs" << endl;
1082 oss << setw(7)<< "" << setw(6) << "Beam"
1083 << setw(23) << dirtype << endl;
1084
1085 oss << asap::SEPERATOR << endl;
1086
1087 // Flush summary and clear up the string
1088 ols << String(oss) << LogIO::POST;
1089 if (ofs) ofs << String(oss) << flush;
1090 oss.str("");
1091 oss.clear();
1092
1093
1094 // Get Freq_ID map
1095 ROScalarColumn<uInt> ftabIds(frequencies().table(), "ID");
1096 Int nfid = ftabIds.nrow();
1097 if (nfid <= 0){
1098 oss << "FREQUENCIES subtable is empty: there are no data!!!" << endl;
1099 oss << asap::SEPERATOR << endl;
1100
1101 ols << String(oss) << LogIO::POST;
1102 if (ofs) {
1103 ofs << String(oss) << flush;
1104 ofs.close();
1105 }
1106 return;
1107 }
1108 // Storages of overall IFNO, POLNO, and nchan per FREQ_ID
1109 // the orders are identical to ID in FREQ subtable
1110 Block< Vector<uInt> > ifNos(nfid), polNos(nfid);
1111 Vector<Int> fIdchans(nfid,-1);
1112 map<uInt, Int> fidMap; // (FREQ_ID, row # in FREQ subtable) pair
1113 for (Int i=0; i < nfid; i++){
1114 // fidMap[freqId] returns row number in FREQ subtable
1115 fidMap.insert(pair<uInt, Int>(ftabIds(i),i));
1116 ifNos[i] = Vector<uInt>();
1117 polNos[i] = Vector<uInt>();
1118 }
1119
1120 TableIterator iter(table_, "SCANNO");
1121
1122 // Vars for keeping track of time, freqids, molIds in a SCANNO
1123 //Vector<uInt> freqids;
1124 //Vector<uInt> molids;
1125 Vector<uInt> beamids(1,0);
1126 Vector<MDirection> beamDirs;
1127 Vector<Int> stypeids(1,0);
1128 Vector<String> stypestrs;
1129 Int nfreq(1);
1130 Int nmol(1);
1131 uInt nbeam(1);
1132 uInt nstype(1);
1133
1134 Double btime(0.0), etime(0.0);
1135 Double meanIntTim(0.0);
1136
1137 uInt currFreqId(0), ftabRow(0);
1138 Int iflen(0), pollen(0);
1139
1140 while (!iter.pastEnd()) {
1141 Table subt = iter.table();
1142 uInt snrow = subt.nrow();
1143 ROTableRow row(subt);
1144 const TableRecord& rec = row.get(0);
1145
1146 // relevant columns
1147 ROScalarColumn<Double> mjdCol(subt,"TIME");
1148 ROScalarColumn<Double> intervalCol(subt,"INTERVAL");
1149 MDirection::ROScalarColumn dirCol(subt,"DIRECTION");
1150
1151 ScalarColumn<uInt> freqIdCol(subt,"FREQ_ID");
1152 ScalarColumn<uInt> molIdCol(subt,"MOLECULE_ID");
1153 ROScalarColumn<uInt> beamCol(subt,"BEAMNO");
1154 ROScalarColumn<Int> stypeCol(subt,"SRCTYPE");
1155
1156 ROScalarColumn<uInt> ifNoCol(subt,"IFNO");
1157 ROScalarColumn<uInt> polNoCol(subt,"POLNO");
1158
1159
1160 // Times
1161 meanIntTim = sum(intervalCol.getColumn()) / (double) snrow;
1162 minMax(btime, etime, mjdCol.getColumn());
1163 etime += meanIntTim/C::day;
1164
1165 // MOLECULE_ID and FREQ_ID
1166 Vector<uInt> molids(getNumbers(molIdCol));
1167 molids.shape(nmol);
1168
1169 Vector<uInt> freqids(getNumbers(freqIdCol));
1170 freqids.shape(nfreq);
1171
1172 // Add first beamid, and srcNames
1173 beamids.resize(1,False);
1174 beamDirs.resize(1,False);
1175 beamids(0)=beamCol(0);
1176 beamDirs(0)=dirCol(0);
1177 nbeam = 1;
1178
1179 stypeids.resize(1,False);
1180 stypeids(0)=stypeCol(0);
1181 nstype = 1;
1182
1183 // Global listings of nchan/IFNO/POLNO per FREQ_ID
1184 currFreqId=freqIdCol(0);
1185 ftabRow = fidMap[currFreqId];
1186 // Assumes an identical number of channels per FREQ_ID
1187 if (fIdchans(ftabRow) < 0 ) {
1188 RORecordFieldPtr< Array<Float> > spec(rec, "SPECTRA");
1189 fIdchans(ftabRow)=(*spec).shape()(0);
1190 }
1191 // Should keep ifNos and polNos form the previous SCANNO
1192 if ( !anyEQ(ifNos[ftabRow],ifNoCol(0)) ) {
1193 ifNos[ftabRow].shape(iflen);
1194 iflen++;
1195 ifNos[ftabRow].resize(iflen,True);
1196 ifNos[ftabRow](iflen-1) = ifNoCol(0);
1197 }
1198 if ( !anyEQ(polNos[ftabRow],polNoCol(0)) ) {
1199 polNos[ftabRow].shape(pollen);
1200 pollen++;
1201 polNos[ftabRow].resize(pollen,True);
1202 polNos[ftabRow](pollen-1) = polNoCol(0);
1203 }
1204
1205 for (uInt i=1; i < snrow; i++){
1206 // Need to list BEAMNO and DIRECTION in the same order
1207 if ( !anyEQ(beamids,beamCol(i)) ) {
1208 nbeam++;
1209 beamids.resize(nbeam,True);
1210 beamids(nbeam-1)=beamCol(i);
1211 beamDirs.resize(nbeam,True);
1212 beamDirs(nbeam-1)=dirCol(i);
1213 }
1214
1215 // SRCTYPE is Int (getNumber takes only uInt)
1216 if ( !anyEQ(stypeids,stypeCol(i)) ) {
1217 nstype++;
1218 stypeids.resize(nstype,True);
1219 stypeids(nstype-1)=stypeCol(i);
1220 }
1221
1222 // Global listings of nchan/IFNO/POLNO per FREQ_ID
1223 currFreqId=freqIdCol(i);
1224 ftabRow = fidMap[currFreqId];
1225 if (fIdchans(ftabRow) < 0 ) {
1226 const TableRecord& rec = row.get(i);
1227 RORecordFieldPtr< Array<Float> > spec(rec, "SPECTRA");
1228 fIdchans(ftabRow) = (*spec).shape()(0);
1229 }
1230 if ( !anyEQ(ifNos[ftabRow],ifNoCol(i)) ) {
1231 ifNos[ftabRow].shape(iflen);
1232 iflen++;
1233 ifNos[ftabRow].resize(iflen,True);
1234 ifNos[ftabRow](iflen-1) = ifNoCol(i);
1235 }
1236 if ( !anyEQ(polNos[ftabRow],polNoCol(i)) ) {
1237 polNos[ftabRow].shape(pollen);
1238 pollen++;
1239 polNos[ftabRow].resize(pollen,True);
1240 polNos[ftabRow](pollen-1) = polNoCol(i);
1241 }
1242 } // end of row iteration
1243
1244 stypestrs.resize(nstype,False);
1245 for (uInt j=0; j < nstype; j++)
1246 stypestrs(j) = SrcType::getName(stypeids(j));
1247
1248 // Format Scan summary
1249 oss << setw(4) << std::right << rec.asuInt("SCANNO")
1250 << std::left << setw(1) << ""
1251 << setw(15) << rec.asString("SRCNAME")
1252 << setw(21) << MVTime(btime).string(MVTime::YMD,7)
1253 << setw(3) << " - " << MVTime(etime).string(MVTime::TIME,7)
1254 << setw(3) << "" << setw(6) << meanIntTim << setw(1) << ""
1255 << std::right << setw(5) << snrow << setw(2) << ""
1256 << std::left << stypestrs << setw(1) << ""
1257 << freqids << setw(1) << ""
1258 << molids << endl;
1259 // Format Beam summary
1260 for (uInt j=0; j < nbeam; j++) {
1261 oss << setw(7) << "" << setw(6) << beamids(j) << setw(1) << ""
1262 << formatDirection(beamDirs(j)) << endl;
1263 }
1264 // Flush summary every scan and clear up the string
1265 ols << String(oss) << LogIO::POST;
1266 if (ofs) ofs << String(oss) << flush;
1267 oss.str("");
1268 oss.clear();
1269
1270 ++iter;
1271 } // end of scan iteration
1272 oss << asap::SEPERATOR << endl;
1273
1274 // List FRECUENCIES Table (using STFrequencies.print may be slow)
1275 oss << "FREQUENCIES: " << nfreq << endl;
1276 oss << std::right << setw(5) << "ID" << setw(2) << ""
1277 << std::left << setw(5) << "IFNO" << setw(2) << ""
1278 << setw(8) << "Frame"
1279 << setw(16) << "RefVal"
1280 << setw(7) << "RefPix"
1281 << setw(15) << "Increment"
1282 << setw(9) << "Channels"
1283 << setw(6) << "POLNOs" << endl;
1284 Int tmplen;
1285 for (Int i=0; i < nfid; i++){
1286 // List row=i of FREQUENCIES subtable
1287 ifNos[i].shape(tmplen);
1288 if (tmplen >= 1) {
1289 oss << std::right << setw(5) << ftabIds(i) << setw(2) << ""
1290 << setw(3) << ifNos[i](0) << setw(1) << ""
1291 << std::left << setw(46) << frequencies().print(ftabIds(i))
1292 << setw(2) << ""
1293 << std::right << setw(8) << fIdchans[i] << setw(2) << ""
1294 << std::left << polNos[i];
1295 if (tmplen > 1) {
1296 oss << " (" << tmplen << " chains)";
1297 }
1298 oss << endl;
1299 }
1300
1301 }
1302 oss << asap::SEPERATOR << endl;
1303
1304 // List MOLECULES Table (currently lists all rows)
1305 oss << "MOLECULES: " << endl;
1306 if (molecules().nrow() <= 0) {
1307 oss << " MOLECULES subtable is empty: there are no data" << endl;
1308 } else {
1309 ROTableRow row(molecules().table());
1310 oss << std::right << setw(5) << "ID"
1311 << std::left << setw(3) << ""
1312 << setw(18) << "RestFreq"
1313 << setw(15) << "Name" << endl;
1314 for (Int i=0; i < molecules().nrow(); i++){
1315 const TableRecord& rec=row.get(i);
1316 oss << std::right << setw(5) << rec.asuInt("ID")
1317 << std::left << setw(3) << ""
1318 << rec.asArrayDouble("RESTFREQUENCY") << setw(1) << ""
1319 << rec.asArrayString("NAME") << endl;
1320 }
1321 }
1322 oss << asap::SEPERATOR << endl;
1323 ols << String(oss) << LogIO::POST;
1324 if (ofs) {
1325 ofs << String(oss) << flush;
1326 ofs.close();
1327 }
1328 // return String(oss);
1329}
1330
1331
1332std::string Scantable::oldheaderSummary()
1333{
1334 // Format header info
1335// STHeader sdh;
1336// sdh = getHeader();
1337// sdh.print();
1338 ostringstream oss;
1339 oss.flags(std::ios_base::left);
1340 oss << setw(15) << "Beams:" << setw(4) << nbeam() << endl
1341 << setw(15) << "IFs:" << setw(4) << nif() << endl
1342 << setw(15) << "Polarisations:" << setw(4) << npol()
1343 << "(" << getPolType() << ")" << endl
1344 << setw(15) << "Channels:" << nchan() << endl;
1345 String tmp;
1346 oss << setw(15) << "Observer:"
1347 << table_.keywordSet().asString("Observer") << endl;
1348 oss << setw(15) << "Obs Date:" << getTime(-1,true) << endl;
1349 table_.keywordSet().get("Project", tmp);
1350 oss << setw(15) << "Project:" << tmp << endl;
1351 table_.keywordSet().get("Obstype", tmp);
1352 oss << setw(15) << "Obs. Type:" << tmp << endl;
1353 table_.keywordSet().get("AntennaName", tmp);
1354 oss << setw(15) << "Antenna Name:" << tmp << endl;
1355 table_.keywordSet().get("FluxUnit", tmp);
1356 oss << setw(15) << "Flux Unit:" << tmp << endl;
1357 int nid = moleculeTable_.nrow();
1358 Bool firstline = True;
1359 oss << setw(15) << "Rest Freqs:";
1360 for (int i=0; i<nid; i++) {
1361 Table t = table_(table_.col("MOLECULE_ID") == i, 1);
1362 if (t.nrow() > 0) {
1363 Vector<Double> vec(moleculeTable_.getRestFrequency(i));
1364 if (vec.nelements() > 0) {
1365 if (firstline) {
1366 oss << setprecision(10) << vec << " [Hz]" << endl;
1367 firstline=False;
1368 }
1369 else{
1370 oss << setw(15)<<" " << setprecision(10) << vec << " [Hz]" << endl;
1371 }
1372 } else {
1373 oss << "none" << endl;
1374 }
1375 }
1376 }
1377
1378 oss << setw(15) << "Abcissa:" << getAbcissaLabel(0) << endl;
1379 oss << selector_.print() << endl;
1380 return String(oss);
1381}
1382
1383 //std::string Scantable::summary( const std::string& filename )
1384void Scantable::oldsummary( const std::string& filename )
1385{
1386 ostringstream oss;
1387 ofstream ofs;
1388 LogIO ols(LogOrigin("Scantable", "summary", WHERE));
1389
1390 if (filename != "")
1391 ofs.open( filename.c_str(), ios::out );
1392
1393 oss << endl;
1394 oss << asap::SEPERATOR << endl;
1395 oss << " Scan Table Summary" << endl;
1396 oss << asap::SEPERATOR << endl;
1397
1398 // Format header info
1399 oss << oldheaderSummary();
1400 oss << endl;
1401
1402 // main table
1403 String dirtype = "Position ("
1404 + getDirectionRefString()
1405 + ")";
1406 oss.flags(std::ios_base::left);
1407 oss << setw(5) << "Scan" << setw(15) << "Source"
1408 << setw(10) << "Time" << setw(18) << "Integration"
1409 << setw(15) << "Source Type" << endl;
1410 oss << setw(5) << "" << setw(5) << "Beam" << setw(3) << "" << dirtype << endl;
1411 oss << setw(10) << "" << setw(3) << "IF" << setw(3) << ""
1412 << setw(8) << "Frame" << setw(16)
1413 << "RefVal" << setw(10) << "RefPix" << setw(12) << "Increment"
1414 << setw(7) << "Channels"
1415 << endl;
1416 oss << asap::SEPERATOR << endl;
1417
1418 // Flush summary and clear up the string
1419 ols << String(oss) << LogIO::POST;
1420 if (ofs) ofs << String(oss) << flush;
1421 oss.str("");
1422 oss.clear();
1423
1424 TableIterator iter(table_, "SCANNO");
1425 while (!iter.pastEnd()) {
1426 Table subt = iter.table();
1427 ROTableRow row(subt);
1428 MEpoch::ROScalarColumn timeCol(subt,"TIME");
1429 const TableRecord& rec = row.get(0);
1430 oss << setw(4) << std::right << rec.asuInt("SCANNO")
1431 << std::left << setw(1) << ""
1432 << setw(15) << rec.asString("SRCNAME")
1433 << setw(10) << formatTime(timeCol(0), false);
1434 // count the cycles in the scan
1435 TableIterator cyciter(subt, "CYCLENO");
1436 int nint = 0;
1437 while (!cyciter.pastEnd()) {
1438 ++nint;
1439 ++cyciter;
1440 }
1441 oss << setw(3) << std::right << nint << setw(3) << " x " << std::left
1442 << setw(11) << formatSec(rec.asFloat("INTERVAL")) << setw(1) << ""
1443 << setw(15) << SrcType::getName(rec.asInt("SRCTYPE")) << endl;
1444
1445 TableIterator biter(subt, "BEAMNO");
1446 while (!biter.pastEnd()) {
1447 Table bsubt = biter.table();
1448 ROTableRow brow(bsubt);
1449 const TableRecord& brec = brow.get(0);
1450 uInt row0 = bsubt.rowNumbers(table_)[0];
1451 oss << setw(5) << "" << setw(4) << std::right << brec.asuInt("BEAMNO")<< std::left;
1452 oss << setw(4) << "" << formatDirection(getDirection(row0)) << endl;
1453 TableIterator iiter(bsubt, "IFNO");
1454 while (!iiter.pastEnd()) {
1455 Table isubt = iiter.table();
1456 ROTableRow irow(isubt);
1457 const TableRecord& irec = irow.get(0);
1458 oss << setw(9) << "";
1459 oss << setw(3) << std::right << irec.asuInt("IFNO") << std::left
1460 << setw(1) << "" << frequencies().print(irec.asuInt("FREQ_ID"))
1461 << setw(3) << "" << nchan(irec.asuInt("IFNO"))
1462 << endl;
1463
1464 ++iiter;
1465 }
1466 ++biter;
1467 }
1468 // Flush summary every scan and clear up the string
1469 ols << String(oss) << LogIO::POST;
1470 if (ofs) ofs << String(oss) << flush;
1471 oss.str("");
1472 oss.clear();
1473
1474 ++iter;
1475 }
1476 oss << asap::SEPERATOR << endl;
1477 ols << String(oss) << LogIO::POST;
1478 if (ofs) {
1479 ofs << String(oss) << flush;
1480 ofs.close();
1481 }
1482 // return String(oss);
1483}
1484
1485// std::string Scantable::getTime(int whichrow, bool showdate) const
1486// {
1487// MEpoch::ROScalarColumn timeCol(table_, "TIME");
1488// MEpoch me;
1489// if (whichrow > -1) {
1490// me = timeCol(uInt(whichrow));
1491// } else {
1492// Double tm;
1493// table_.keywordSet().get("UTC",tm);
1494// me = MEpoch(MVEpoch(tm));
1495// }
1496// return formatTime(me, showdate);
1497// }
1498
1499std::string Scantable::getTime(int whichrow, bool showdate, uInt prec) const
1500{
1501 MEpoch me;
1502 me = getEpoch(whichrow);
1503 return formatTime(me, showdate, prec);
1504}
1505
1506MEpoch Scantable::getEpoch(int whichrow) const
1507{
1508 if (whichrow > -1) {
1509 return timeCol_(uInt(whichrow));
1510 } else {
1511 Double tm;
1512 table_.keywordSet().get("UTC",tm);
1513 return MEpoch(MVEpoch(tm));
1514 }
1515}
1516
1517std::string Scantable::getDirectionString(int whichrow) const
1518{
1519 return formatDirection(getDirection(uInt(whichrow)));
1520}
1521
1522
1523SpectralCoordinate Scantable::getSpectralCoordinate(int whichrow) const {
1524 const MPosition& mp = getAntennaPosition();
1525 const MDirection& md = getDirection(whichrow);
1526 const MEpoch& me = timeCol_(whichrow);
1527 //Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
1528 Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
1529 return freqTable_.getSpectralCoordinate(md, mp, me, rf,
1530 mfreqidCol_(whichrow));
1531}
1532
1533std::vector< double > Scantable::getAbcissa( int whichrow ) const
1534{
1535 if ( whichrow > int(table_.nrow()) ) throw(AipsError("Illegal row number"));
1536 std::vector<double> stlout;
1537 int nchan = specCol_(whichrow).nelements();
1538 String us = freqTable_.getUnitString();
1539 if ( us == "" || us == "pixel" || us == "channel" ) {
1540 for (int i=0; i<nchan; ++i) {
1541 stlout.push_back(double(i));
1542 }
1543 return stlout;
1544 }
1545 SpectralCoordinate spc = getSpectralCoordinate(whichrow);
1546 Vector<Double> pixel(nchan);
1547 Vector<Double> world;
1548 indgen(pixel);
1549 if ( Unit(us) == Unit("Hz") ) {
1550 for ( int i=0; i < nchan; ++i) {
1551 Double world;
1552 spc.toWorld(world, pixel[i]);
1553 stlout.push_back(double(world));
1554 }
1555 } else if ( Unit(us) == Unit("km/s") ) {
1556 Vector<Double> world;
1557 spc.pixelToVelocity(world, pixel);
1558 world.tovector(stlout);
1559 }
1560 return stlout;
1561}
1562void Scantable::setDirectionRefString( const std::string & refstr )
1563{
1564 MDirection::Types mdt;
1565 if (refstr != "" && !MDirection::getType(mdt, refstr)) {
1566 throw(AipsError("Illegal Direction frame."));
1567 }
1568 if ( refstr == "" ) {
1569 String defaultstr = MDirection::showType(dirCol_.getMeasRef().getType());
1570 table_.rwKeywordSet().define("DIRECTIONREF", defaultstr);
1571 } else {
1572 table_.rwKeywordSet().define("DIRECTIONREF", String(refstr));
1573 }
1574}
1575
1576std::string Scantable::getDirectionRefString( ) const
1577{
1578 return table_.keywordSet().asString("DIRECTIONREF");
1579}
1580
1581MDirection Scantable::getDirection(int whichrow ) const
1582{
1583 String usertype = table_.keywordSet().asString("DIRECTIONREF");
1584 String type = MDirection::showType(dirCol_.getMeasRef().getType());
1585 if ( usertype != type ) {
1586 MDirection::Types mdt;
1587 if (!MDirection::getType(mdt, usertype)) {
1588 throw(AipsError("Illegal Direction frame."));
1589 }
1590 return dirCol_.convert(uInt(whichrow), mdt);
1591 } else {
1592 return dirCol_(uInt(whichrow));
1593 }
1594}
1595
1596std::string Scantable::getAbcissaLabel( int whichrow ) const
1597{
1598 if ( whichrow > int(table_.nrow()) ) throw(AipsError("Illegal ro number"));
1599 const MPosition& mp = getAntennaPosition();
1600 const MDirection& md = getDirection(whichrow);
1601 const MEpoch& me = timeCol_(whichrow);
1602 //const Double& rf = mmolidCol_(whichrow);
1603 const Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
1604 SpectralCoordinate spc =
1605 freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
1606
1607 String s = "Channel";
1608 Unit u = Unit(freqTable_.getUnitString());
1609 if (u == Unit("km/s")) {
1610 s = CoordinateUtil::axisLabel(spc, 0, True,True, True);
1611 } else if (u == Unit("Hz")) {
1612 Vector<String> wau(1);wau = u.getName();
1613 spc.setWorldAxisUnits(wau);
1614 s = CoordinateUtil::axisLabel(spc, 0, True, True, False);
1615 }
1616 return s;
1617
1618}
1619
1620/**
1621void asap::Scantable::setRestFrequencies( double rf, const std::string& name,
1622 const std::string& unit )
1623**/
1624void Scantable::setRestFrequencies( vector<double> rf, const vector<std::string>& name,
1625 const std::string& unit )
1626
1627{
1628 ///@todo lookup in line table to fill in name and formattedname
1629 Unit u(unit);
1630 //Quantum<Double> urf(rf, u);
1631 Quantum<Vector<Double> >urf(rf, u);
1632 Vector<String> formattedname(0);
1633 //cerr<<"Scantable::setRestFrequnecies="<<urf<<endl;
1634
1635 //uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), name, "");
1636 uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), mathutil::toVectorString(name), formattedname);
1637 TableVector<uInt> tabvec(table_, "MOLECULE_ID");
1638 tabvec = id;
1639}
1640
1641/**
1642void asap::Scantable::setRestFrequencies( const std::string& name )
1643{
1644 throw(AipsError("setRestFrequencies( const std::string& name ) NYI"));
1645 ///@todo implement
1646}
1647**/
1648
1649void Scantable::setRestFrequencies( const vector<std::string>& name )
1650{
1651 (void) name; // suppress unused warning
1652 throw(AipsError("setRestFrequencies( const vector<std::string>& name ) NYI"));
1653 ///@todo implement
1654}
1655
1656std::vector< unsigned int > Scantable::rownumbers( ) const
1657{
1658 std::vector<unsigned int> stlout;
1659 Vector<uInt> vec = table_.rowNumbers();
1660 vec.tovector(stlout);
1661 return stlout;
1662}
1663
1664
1665Matrix<Float> Scantable::getPolMatrix( uInt whichrow ) const
1666{
1667 ROTableRow row(table_);
1668 const TableRecord& rec = row.get(whichrow);
1669 Table t =
1670 originalTable_( originalTable_.col("SCANNO") == Int(rec.asuInt("SCANNO"))
1671 && originalTable_.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
1672 && originalTable_.col("IFNO") == Int(rec.asuInt("IFNO"))
1673 && originalTable_.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
1674 ROArrayColumn<Float> speccol(t, "SPECTRA");
1675 return speccol.getColumn();
1676}
1677
1678std::vector< std::string > Scantable::columnNames( ) const
1679{
1680 Vector<String> vec = table_.tableDesc().columnNames();
1681 return mathutil::tovectorstring(vec);
1682}
1683
1684MEpoch::Types Scantable::getTimeReference( ) const
1685{
1686 return MEpoch::castType(timeCol_.getMeasRef().getType());
1687}
1688
1689void Scantable::addFit( const STFitEntry& fit, int row )
1690{
1691 //cout << mfitidCol_(uInt(row)) << endl;
1692 LogIO os( LogOrigin( "Scantable", "addFit()", WHERE ) ) ;
1693 os << mfitidCol_(uInt(row)) << LogIO::POST ;
1694 uInt id = fitTable_.addEntry(fit, mfitidCol_(uInt(row)));
1695 mfitidCol_.put(uInt(row), id);
1696}
1697
1698void Scantable::shift(int npix)
1699{
1700 Vector<uInt> fids(mfreqidCol_.getColumn());
1701 genSort( fids, Sort::Ascending,
1702 Sort::QuickSort|Sort::NoDuplicates );
1703 for (uInt i=0; i<fids.nelements(); ++i) {
1704 frequencies().shiftRefPix(npix, fids[i]);
1705 }
1706}
1707
1708String Scantable::getAntennaName() const
1709{
1710 String out;
1711 table_.keywordSet().get("AntennaName", out);
1712 String::size_type pos1 = out.find("@") ;
1713 String::size_type pos2 = out.find("//") ;
1714 if ( pos2 != String::npos )
1715 out = out.substr(pos2+2,pos1-pos2-2) ;
1716 else if ( pos1 != String::npos )
1717 out = out.substr(0,pos1) ;
1718 return out;
1719}
1720
1721int Scantable::checkScanInfo(const std::vector<int>& scanlist) const
1722{
1723 String tbpath;
1724 int ret = 0;
1725 if ( table_.keywordSet().isDefined("GBT_GO") ) {
1726 table_.keywordSet().get("GBT_GO", tbpath);
1727 Table t(tbpath,Table::Old);
1728 // check each scan if other scan of the pair exist
1729 int nscan = scanlist.size();
1730 for (int i = 0; i < nscan; i++) {
1731 Table subt = t( t.col("SCAN") == scanlist[i] );
1732 if (subt.nrow()==0) {
1733 //cerr <<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<endl;
1734 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1735 os <<LogIO::WARN<<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<LogIO::POST;
1736 ret = 1;
1737 break;
1738 }
1739 ROTableRow row(subt);
1740 const TableRecord& rec = row.get(0);
1741 int scan1seqn = rec.asuInt("PROCSEQN");
1742 int laston1 = rec.asuInt("LASTON");
1743 if ( rec.asuInt("PROCSIZE")==2 ) {
1744 if ( i < nscan-1 ) {
1745 Table subt2 = t( t.col("SCAN") == scanlist[i+1] );
1746 if ( subt2.nrow() == 0) {
1747 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1748
1749 //cerr<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<endl;
1750 os<<LogIO::WARN<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<LogIO::POST;
1751 ret = 1;
1752 break;
1753 }
1754 ROTableRow row2(subt2);
1755 const TableRecord& rec2 = row2.get(0);
1756 int scan2seqn = rec2.asuInt("PROCSEQN");
1757 int laston2 = rec2.asuInt("LASTON");
1758 if (scan1seqn == 1 && scan2seqn == 2) {
1759 if (laston1 == laston2) {
1760 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1761 //cerr<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
1762 os<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST;
1763 i +=1;
1764 }
1765 else {
1766 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1767 //cerr<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
1768 os<<LogIO::WARN<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST;
1769 }
1770 }
1771 else if (scan1seqn==2 && scan2seqn == 1) {
1772 if (laston1 == laston2) {
1773 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1774 //cerr<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<endl;
1775 os<<LogIO::WARN<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<LogIO::POST;
1776 ret = 1;
1777 break;
1778 }
1779 }
1780 else {
1781 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1782 //cerr<<"The other scan for "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<endl;
1783 os<<LogIO::WARN<<"The other scan for "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<LogIO::POST;
1784 ret = 1;
1785 break;
1786 }
1787 }
1788 }
1789 else {
1790 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1791 //cerr<<"The scan does not appear to be standard obsevation."<<endl;
1792 os<<LogIO::WARN<<"The scan does not appear to be standard obsevation."<<LogIO::POST;
1793 }
1794 //if ( i >= nscan ) break;
1795 }
1796 }
1797 else {
1798 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1799 //cerr<<"No reference to GBT_GO table."<<endl;
1800 os<<LogIO::WARN<<"No reference to GBT_GO table."<<LogIO::POST;
1801 ret = 1;
1802 }
1803 return ret;
1804}
1805
1806std::vector<double> Scantable::getDirectionVector(int whichrow) const
1807{
1808 Vector<Double> Dir = dirCol_(whichrow).getAngle("rad").getValue();
1809 std::vector<double> dir;
1810 Dir.tovector(dir);
1811 return dir;
1812}
1813
1814void asap::Scantable::reshapeSpectrum( int nmin, int nmax )
1815 throw( casa::AipsError )
1816{
1817 // assumed that all rows have same nChan
1818 Vector<Float> arr = specCol_( 0 ) ;
1819 int nChan = arr.nelements() ;
1820
1821 // if nmin < 0 or nmax < 0, nothing to do
1822 if ( nmin < 0 ) {
1823 throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
1824 }
1825 if ( nmax < 0 ) {
1826 throw( casa::indexError<int>( nmax, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
1827 }
1828
1829 // if nmin > nmax, exchange values
1830 if ( nmin > nmax ) {
1831 int tmp = nmax ;
1832 nmax = nmin ;
1833 nmin = tmp ;
1834 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1835 os << "Swap values. Applied range is ["
1836 << nmin << ", " << nmax << "]" << LogIO::POST ;
1837 }
1838
1839 // if nmin exceeds nChan, nothing to do
1840 if ( nmin >= nChan ) {
1841 throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Specified minimum exceeds nChan." ) ) ;
1842 }
1843
1844 // if nmax exceeds nChan, reset nmax to nChan
1845 if ( nmax >= nChan-1 ) {
1846 if ( nmin == 0 ) {
1847 // nothing to do
1848 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1849 os << "Whole range is selected. Nothing to do." << LogIO::POST ;
1850 return ;
1851 }
1852 else {
1853 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1854 os << "Specified maximum exceeds nChan. Applied range is ["
1855 << nmin << ", " << nChan-1 << "]." << LogIO::POST ;
1856 nmax = nChan - 1 ;
1857 }
1858 }
1859
1860 // reshape specCol_ and flagCol_
1861 for ( int irow = 0 ; irow < nrow() ; irow++ ) {
1862 reshapeSpectrum( nmin, nmax, irow ) ;
1863 }
1864
1865 // update FREQUENCIES subtable
1866 Double refpix ;
1867 Double refval ;
1868 Double increment ;
1869 int freqnrow = freqTable_.table().nrow() ;
1870 Vector<uInt> oldId( freqnrow ) ;
1871 Vector<uInt> newId( freqnrow ) ;
1872 for ( int irow = 0 ; irow < freqnrow ; irow++ ) {
1873 freqTable_.getEntry( refpix, refval, increment, irow ) ;
1874 /***
1875 * need to shift refpix to nmin
1876 * note that channel nmin in old index will be channel 0 in new one
1877 ***/
1878 refval = refval - ( refpix - nmin ) * increment ;
1879 refpix = 0 ;
1880 freqTable_.setEntry( refpix, refval, increment, irow ) ;
1881 }
1882
1883 // update nchan
1884 int newsize = nmax - nmin + 1 ;
1885 table_.rwKeywordSet().define( "nChan", newsize ) ;
1886
1887 // update bandwidth
1888 // assumed all spectra in the scantable have same bandwidth
1889 table_.rwKeywordSet().define( "Bandwidth", increment * newsize ) ;
1890
1891 return ;
1892}
1893
1894void asap::Scantable::reshapeSpectrum( int nmin, int nmax, int irow )
1895{
1896 // reshape specCol_ and flagCol_
1897 Vector<Float> oldspec = specCol_( irow ) ;
1898 Vector<uChar> oldflag = flagsCol_( irow ) ;
1899 Vector<Float> oldtsys = tsysCol_( irow ) ;
1900 uInt newsize = nmax - nmin + 1 ;
1901 Slice slice( nmin, newsize, 1 ) ;
1902 specCol_.put( irow, oldspec( slice ) ) ;
1903 flagsCol_.put( irow, oldflag( slice ) ) ;
1904 if ( oldspec.size() == oldtsys.size() )
1905 tsysCol_.put( irow, oldtsys( slice ) ) ;
1906
1907 return ;
1908}
1909
1910void asap::Scantable::regridSpecChannel( double dnu, int nChan )
1911{
1912 LogIO os( LogOrigin( "Scantable", "regridChannel()", WHERE ) ) ;
1913 os << "Regrid abcissa with spectral resoultion " << dnu << " " << freqTable_.getUnitString() << " with channel number " << ((nChan>0)? String(nChan) : "covering band width")<< LogIO::POST ;
1914 int freqnrow = freqTable_.table().nrow() ;
1915 Vector<bool> firstTime( freqnrow, true ) ;
1916 double oldincr, factor;
1917 uInt currId;
1918 Double refpix ;
1919 Double refval ;
1920 Double increment ;
1921 for ( int irow = 0 ; irow < nrow() ; irow++ ) {
1922 currId = mfreqidCol_(irow);
1923 vector<double> abcissa = getAbcissa( irow ) ;
1924 if (nChan < 0) {
1925 int oldsize = abcissa.size() ;
1926 double bw = (abcissa[oldsize-1]-abcissa[0]) + \
1927 0.5 * (abcissa[1]-abcissa[0] + abcissa[oldsize-1]-abcissa[oldsize-2]) ;
1928 nChan = int( ceil( abs(bw/dnu) ) ) ;
1929 }
1930 // actual regridding
1931 regridChannel( nChan, dnu, irow ) ;
1932
1933 // update FREQUENCIES subtable
1934 if (firstTime[currId]) {
1935 oldincr = abcissa[1]-abcissa[0] ;
1936 factor = dnu/oldincr ;
1937 firstTime[currId] = false ;
1938 freqTable_.getEntry( refpix, refval, increment, currId ) ;
1939
1940 //refval = refval - ( refpix + 0.5 * (1 - factor) ) * increment ;
1941 if (factor > 0 ) {
1942 refpix = (refpix + 0.5)/factor - 0.5;
1943 } else {
1944 refpix = (abcissa.size() - 0.5 - refpix)/abs(factor) - 0.5;
1945 }
1946 freqTable_.setEntry( refpix, refval, increment*factor, currId ) ;
1947 //os << "ID" << currId << ": channel width (Orig) = " << oldincr << " [" << freqTable_.getUnitString() << "], scale factor = " << factor << LogIO::POST ;
1948 //os << " frequency increment (Orig) = " << increment << "-> (New) " << increment*factor << LogIO::POST ;
1949 }
1950 }
1951}
1952
1953void asap::Scantable::regridChannel( int nChan, double dnu )
1954{
1955 LogIO os( LogOrigin( "Scantable", "regridChannel()", WHERE ) ) ;
1956 os << "Regrid abcissa with channel number " << nChan << " and spectral resoultion " << dnu << "Hz." << LogIO::POST ;
1957 // assumed that all rows have same nChan
1958 Vector<Float> arr = specCol_( 0 ) ;
1959 int oldsize = arr.nelements() ;
1960
1961 // if oldsize == nChan, nothing to do
1962 if ( oldsize == nChan ) {
1963 os << "Specified channel number is same as current one. Nothing to do." << LogIO::POST ;
1964 return ;
1965 }
1966
1967 // if oldChan < nChan, unphysical operation
1968 if ( oldsize < nChan ) {
1969 os << "Unphysical operation. Nothing to do." << LogIO::POST ;
1970 return ;
1971 }
1972
1973 // change channel number for specCol_, flagCol_, and tsysCol_ (if necessary)
1974 vector<string> coordinfo = getCoordInfo() ;
1975 string oldinfo = coordinfo[0] ;
1976 coordinfo[0] = "Hz" ;
1977 setCoordInfo( coordinfo ) ;
1978 for ( int irow = 0 ; irow < nrow() ; irow++ ) {
1979 regridChannel( nChan, dnu, irow ) ;
1980 }
1981 coordinfo[0] = oldinfo ;
1982 setCoordInfo( coordinfo ) ;
1983
1984
1985 // NOTE: this method does not update metadata such as
1986 // FREQUENCIES subtable, nChan, Bandwidth, etc.
1987
1988 return ;
1989}
1990
1991void asap::Scantable::regridChannel( int nChan, double dnu, int irow )
1992{
1993 // logging
1994 //ofstream ofs( "average.log", std::ios::out | std::ios::app ) ;
1995 //ofs << "IFNO = " << getIF( irow ) << " irow = " << irow << endl ;
1996
1997 Vector<Float> oldspec = specCol_( irow ) ;
1998 Vector<uChar> oldflag = flagsCol_( irow ) ;
1999 Vector<Float> oldtsys = tsysCol_( irow ) ;
2000 Vector<Float> newspec( nChan, 0 ) ;
2001 Vector<uChar> newflag( nChan, true ) ;
2002 Vector<Float> newtsys ;
2003 bool regridTsys = false ;
2004 if (oldtsys.size() == oldspec.size()) {
2005 regridTsys = true ;
2006 newtsys.resize(nChan,false) ;
2007 newtsys = 0 ;
2008 }
2009
2010 // regrid
2011 vector<double> abcissa = getAbcissa( irow ) ;
2012 int oldsize = abcissa.size() ;
2013 double olddnu = abcissa[1] - abcissa[0] ;
2014 //int ichan = 0 ;
2015 double wsum = 0.0 ;
2016 Vector<double> zi( nChan+1 ) ;
2017 Vector<double> yi( oldsize + 1 ) ;
2018 yi[0] = abcissa[0] - 0.5 * olddnu ;
2019 for ( int ii = 1 ; ii < oldsize ; ii++ )
2020 yi[ii] = 0.5* (abcissa[ii-1] + abcissa[ii]) ;
2021 yi[oldsize] = abcissa[oldsize-1] \
2022 + 0.5 * (abcissa[oldsize-1] - abcissa[oldsize-2]) ;
2023 //zi[0] = abcissa[0] - 0.5 * olddnu ;
2024 zi[0] = ((olddnu*dnu > 0) ? yi[0] : yi[oldsize]) ;
2025 for ( int ii = 1 ; ii < nChan ; ii++ )
2026 zi[ii] = zi[0] + dnu * ii ;
2027 zi[nChan] = zi[nChan-1] + dnu ;
2028 // Access zi and yi in ascending order
2029 int izs = ((dnu > 0) ? 0 : nChan ) ;
2030 int ize = ((dnu > 0) ? nChan : 0 ) ;
2031 int izincr = ((dnu > 0) ? 1 : -1 ) ;
2032 int ichan = ((olddnu > 0) ? 0 : oldsize ) ;
2033 int iye = ((olddnu > 0) ? oldsize : 0 ) ;
2034 int iyincr = ((olddnu > 0) ? 1 : -1 ) ;
2035 //for ( int ii = izs ; ii != ize ; ii+=izincr ){
2036 int ii = izs ;
2037 while (ii != ize) {
2038 // always zl < zr
2039 double zl = zi[ii] ;
2040 double zr = zi[ii+izincr] ;
2041 // Need to access smaller index for the new spec, flag, and tsys.
2042 // Values between zi[k] and zi[k+1] should be stored in newspec[k], etc.
2043 int i = min(ii, ii+izincr) ;
2044 //for ( int jj = ichan ; jj != iye ; jj+=iyincr ) {
2045 int jj = ichan ;
2046 while (jj != iye) {
2047 // always yl < yr
2048 double yl = yi[jj] ;
2049 double yr = yi[jj+iyincr] ;
2050 // Need to access smaller index for the original spec, flag, and tsys.
2051 // Values between yi[k] and yi[k+1] are stored in oldspec[k], etc.
2052 int j = min(jj, jj+iyincr) ;
2053 if ( yr <= zl ) {
2054 jj += iyincr ;
2055 continue ;
2056 }
2057 else if ( yl <= zl ) {
2058 if ( yr < zr ) {
2059 if (!oldflag[j]) {
2060 newspec[i] += oldspec[j] * ( yr - zl ) ;
2061 if (regridTsys) newtsys[i] += oldtsys[j] * ( yr - zl ) ;
2062 wsum += ( yr - zl ) ;
2063 }
2064 newflag[i] = newflag[i] && oldflag[j] ;
2065 }
2066 else {
2067 if (!oldflag[j]) {
2068 newspec[i] += oldspec[j] * abs(dnu) ;
2069 if (regridTsys) newtsys[i] += oldtsys[j] * abs(dnu) ;
2070 wsum += abs(dnu) ;
2071 }
2072 newflag[i] = newflag[i] && oldflag[j] ;
2073 ichan = jj ;
2074 break ;
2075 }
2076 }
2077 else if ( yl < zr ) {
2078 if ( yr <= zr ) {
2079 if (!oldflag[j]) {
2080 newspec[i] += oldspec[j] * ( yr - yl ) ;
2081 if (regridTsys) newtsys[i] += oldtsys[j] * ( yr - yl ) ;
2082 wsum += ( yr - yl ) ;
2083 }
2084 newflag[i] = newflag[i] && oldflag[j] ;
2085 }
2086 else {
2087 if (!oldflag[j]) {
2088 newspec[i] += oldspec[j] * ( zr - yl ) ;
2089 if (regridTsys) newtsys[i] += oldtsys[j] * ( zr - yl ) ;
2090 wsum += ( zr - yl ) ;
2091 }
2092 newflag[i] = newflag[i] && oldflag[j] ;
2093 ichan = jj ;
2094 break ;
2095 }
2096 }
2097 else {
2098 ichan = jj - iyincr ;
2099 break ;
2100 }
2101 jj += iyincr ;
2102 }
2103 if ( wsum != 0.0 ) {
2104 newspec[i] /= wsum ;
2105 if (regridTsys) newtsys[i] /= wsum ;
2106 }
2107 wsum = 0.0 ;
2108 ii += izincr ;
2109 }
2110// if ( dnu > 0.0 ) {
2111// for ( int ii = 0 ; ii < nChan ; ii++ ) {
2112// double zl = zi[ii] ;
2113// double zr = zi[ii+1] ;
2114// for ( int j = ichan ; j < oldsize ; j++ ) {
2115// double yl = yi[j] ;
2116// double yr = yi[j+1] ;
2117// if ( yl <= zl ) {
2118// if ( yr <= zl ) {
2119// continue ;
2120// }
2121// else if ( yr <= zr ) {
2122// if (!oldflag[j]) {
2123// newspec[ii] += oldspec[j] * ( yr - zl ) ;
2124// if (regridTsys) newtsys[ii] += oldtsys[j] * ( yr - zl ) ;
2125// wsum += ( yr - zl ) ;
2126// }
2127// newflag[ii] = newflag[ii] && oldflag[j] ;
2128// }
2129// else {
2130// if (!oldflag[j]) {
2131// newspec[ii] += oldspec[j] * dnu ;
2132// if (regridTsys) newtsys[ii] += oldtsys[j] * dnu ;
2133// wsum += dnu ;
2134// }
2135// newflag[ii] = newflag[ii] && oldflag[j] ;
2136// ichan = j ;
2137// break ;
2138// }
2139// }
2140// else if ( yl < zr ) {
2141// if ( yr <= zr ) {
2142// if (!oldflag[j]) {
2143// newspec[ii] += oldspec[j] * ( yr - yl ) ;
2144// if (regridTsys) newtsys[ii] += oldtsys[j] * ( yr - yl ) ;
2145// wsum += ( yr - yl ) ;
2146// }
2147// newflag[ii] = newflag[ii] && oldflag[j] ;
2148// }
2149// else {
2150// if (!oldflag[j]) {
2151// newspec[ii] += oldspec[j] * ( zr - yl ) ;
2152// if (regridTsys) newtsys[ii] += oldtsys[j] * ( zr - yl ) ;
2153// wsum += ( zr - yl ) ;
2154// }
2155// newflag[ii] = newflag[ii] && oldflag[j] ;
2156// ichan = j ;
2157// break ;
2158// }
2159// }
2160// else {
2161// ichan = j - 1 ;
2162// break ;
2163// }
2164// }
2165// if ( wsum != 0.0 ) {
2166// newspec[ii] /= wsum ;
2167// if (regridTsys) newtsys[ii] /= wsum ;
2168// }
2169// wsum = 0.0 ;
2170// }
2171// }
2172// else if ( dnu < 0.0 ) {
2173// for ( int ii = 0 ; ii < nChan ; ii++ ) {
2174// double zl = zi[ii] ;
2175// double zr = zi[ii+1] ;
2176// for ( int j = ichan ; j < oldsize ; j++ ) {
2177// double yl = yi[j] ;
2178// double yr = yi[j+1] ;
2179// if ( yl >= zl ) {
2180// if ( yr >= zl ) {
2181// continue ;
2182// }
2183// else if ( yr >= zr ) {
2184// if (!oldflag[j]) {
2185// newspec[ii] += oldspec[j] * abs( yr - zl ) ;
2186// if (regridTsys) newtsys[ii] += oldtsys[j] * abs( yr - zl ) ;
2187// wsum += abs( yr - zl ) ;
2188// }
2189// newflag[ii] = newflag[ii] && oldflag[j] ;
2190// }
2191// else {
2192// if (!oldflag[j]) {
2193// newspec[ii] += oldspec[j] * abs( dnu ) ;
2194// if (regridTsys) newtsys[ii] += oldtsys[j] * abs( dnu ) ;
2195// wsum += abs( dnu ) ;
2196// }
2197// newflag[ii] = newflag[ii] && oldflag[j] ;
2198// ichan = j ;
2199// break ;
2200// }
2201// }
2202// else if ( yl > zr ) {
2203// if ( yr >= zr ) {
2204// if (!oldflag[j]) {
2205// newspec[ii] += oldspec[j] * abs( yr - yl ) ;
2206// if (regridTsys) newtsys[ii] += oldtsys[j] * abs( yr - yl ) ;
2207// wsum += abs( yr - yl ) ;
2208// }
2209// newflag[ii] = newflag[ii] && oldflag[j] ;
2210// }
2211// else {
2212// if (!oldflag[j]) {
2213// newspec[ii] += oldspec[j] * abs( zr - yl ) ;
2214// if (regridTsys) newtsys[ii] += oldtsys[j] * abs( zr - yl ) ;
2215// wsum += abs( zr - yl ) ;
2216// }
2217// newflag[ii] = newflag[ii] && oldflag[j] ;
2218// ichan = j ;
2219// break ;
2220// }
2221// }
2222// else {
2223// ichan = j - 1 ;
2224// break ;
2225// }
2226// }
2227// if ( wsum != 0.0 ) {
2228// newspec[ii] /= wsum ;
2229// if (regridTsys) newtsys[ii] /= wsum ;
2230// }
2231// wsum = 0.0 ;
2232// }
2233// }
2234// // //ofs << "olddnu = " << olddnu << ", dnu = " << dnu << endl ;
2235// // pile += dnu ;
2236// // wedge = olddnu * ( refChan + 1 ) ;
2237// // while ( wedge < pile ) {
2238// // newspec[0] += olddnu * oldspec[refChan] ;
2239// // newflag[0] = newflag[0] || oldflag[refChan] ;
2240// // //ofs << "channel " << refChan << " is included in new channel 0" << endl ;
2241// // refChan++ ;
2242// // wedge += olddnu ;
2243// // wsum += olddnu ;
2244// // //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
2245// // }
2246// // frac = ( wedge - pile ) / olddnu ;
2247// // wsum += ( 1.0 - frac ) * olddnu ;
2248// // newspec[0] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
2249// // newflag[0] = newflag[0] || oldflag[refChan] ;
2250// // //ofs << "channel " << refChan << " is partly included in new channel 0" << " with fraction of " << ( 1.0 - frac ) << endl ;
2251// // //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
2252// // newspec[0] /= wsum ;
2253// // //ofs << "newspec[0] = " << newspec[0] << endl ;
2254// // //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
2255
2256// // /***
2257// // * ichan = 1 - nChan-2
2258// // ***/
2259// // for ( int ichan = 1 ; ichan < nChan - 1 ; ichan++ ) {
2260// // pile += dnu ;
2261// // newspec[ichan] += frac * olddnu * oldspec[refChan] ;
2262// // newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
2263// // //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << frac << endl ;
2264// // refChan++ ;
2265// // wedge += olddnu ;
2266// // wsum = frac * olddnu ;
2267// // //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
2268// // while ( wedge < pile ) {
2269// // newspec[ichan] += olddnu * oldspec[refChan] ;
2270// // newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
2271// // //ofs << "channel " << refChan << " is included in new channel " << ichan << endl ;
2272// // refChan++ ;
2273// // wedge += olddnu ;
2274// // wsum += olddnu ;
2275// // //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
2276// // }
2277// // frac = ( wedge - pile ) / olddnu ;
2278// // wsum += ( 1.0 - frac ) * olddnu ;
2279// // newspec[ichan] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
2280// // newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
2281// // //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << ( 1.0 - frac ) << endl ;
2282// // //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
2283// // //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
2284// // newspec[ichan] /= wsum ;
2285// // //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << endl ;
2286// // }
2287
2288// // /***
2289// // * ichan = nChan-1
2290// // ***/
2291// // // NOTE: Assumed that all spectra have the same bandwidth
2292// // pile += dnu ;
2293// // newspec[nChan-1] += frac * olddnu * oldspec[refChan] ;
2294// // newflag[nChan-1] = newflag[nChan-1] || oldflag[refChan] ;
2295// // //ofs << "channel " << refChan << " is partly included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
2296// // refChan++ ;
2297// // wedge += olddnu ;
2298// // wsum = frac * olddnu ;
2299// // //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
2300// // for ( int jchan = refChan ; jchan < oldsize ; jchan++ ) {
2301// // newspec[nChan-1] += olddnu * oldspec[jchan] ;
2302// // newflag[nChan-1] = newflag[nChan-1] || oldflag[jchan] ;
2303// // wsum += olddnu ;
2304// // //ofs << "channel " << jchan << " is included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
2305// // //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
2306// // }
2307// // //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
2308// // //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
2309// // newspec[nChan-1] /= wsum ;
2310// // //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << endl ;
2311
2312// // // ofs.close() ;
2313
2314 specCol_.put( irow, newspec ) ;
2315 flagsCol_.put( irow, newflag ) ;
2316 if (regridTsys) tsysCol_.put( irow, newtsys );
2317
2318 return ;
2319}
2320
2321void Scantable::regridChannel( int nChan, double dnu, double fmin, int irow )
2322{
2323 Vector<Float> oldspec = specCol_( irow ) ;
2324 Vector<uChar> oldflag = flagsCol_( irow ) ;
2325 Vector<Float> oldtsys = tsysCol_( irow ) ;
2326 Vector<Float> newspec( nChan, 0 ) ;
2327 Vector<uChar> newflag( nChan, true ) ;
2328 Vector<Float> newtsys ;
2329 bool regridTsys = false ;
2330 if (oldtsys.size() == oldspec.size()) {
2331 regridTsys = true ;
2332 newtsys.resize(nChan,false) ;
2333 newtsys = 0 ;
2334 }
2335
2336 // regrid
2337 vector<double> abcissa = getAbcissa( irow ) ;
2338 int oldsize = abcissa.size() ;
2339 double olddnu = abcissa[1] - abcissa[0] ;
2340 //int ichan = 0 ;
2341 double wsum = 0.0 ;
2342 Vector<double> zi( nChan+1 ) ;
2343 Vector<double> yi( oldsize + 1 ) ;
2344 Block<uInt> count( nChan, 0 ) ;
2345 yi[0] = abcissa[0] - 0.5 * olddnu ;
2346 for ( int ii = 1 ; ii < oldsize ; ii++ )
2347 yi[ii] = 0.5* (abcissa[ii-1] + abcissa[ii]) ;
2348 yi[oldsize] = abcissa[oldsize-1] \
2349 + 0.5 * (abcissa[oldsize-1] - abcissa[oldsize-2]) ;
2350// cout << "olddnu=" << olddnu << ", dnu=" << dnu << " (diff=" << olddnu-dnu << ")" << endl ;
2351// cout << "yi[0]=" << yi[0] << ", fmin=" << fmin << " (diff=" << yi[0]-fmin << ")" << endl ;
2352// cout << "oldsize=" << oldsize << ", nChan=" << nChan << endl ;
2353
2354 // do not regrid if input parameters are almost same as current
2355 // spectral setup
2356 double dnuDiff = abs( ( dnu - olddnu ) / olddnu ) ;
2357 double oldfmin = min( yi[0], yi[oldsize] ) ;
2358 double fminDiff = abs( ( fmin - oldfmin ) / oldfmin ) ;
2359 double nChanDiff = nChan - oldsize ;
2360 double eps = 1.0e-8 ;
2361 if ( nChanDiff == 0 && dnuDiff < eps && fminDiff < eps )
2362 return ;
2363
2364 //zi[0] = abcissa[0] - 0.5 * olddnu ;
2365 //zi[0] = ((olddnu*dnu > 0) ? yi[0] : yi[oldsize]) ;
2366 if ( dnu > 0 )
2367 zi[0] = fmin - 0.5 * dnu ;
2368 else
2369 zi[0] = fmin + nChan * abs(dnu) ;
2370 for ( int ii = 1 ; ii < nChan ; ii++ )
2371 zi[ii] = zi[0] + dnu * ii ;
2372 zi[nChan] = zi[nChan-1] + dnu ;
2373 // Access zi and yi in ascending order
2374 int izs = ((dnu > 0) ? 0 : nChan ) ;
2375 int ize = ((dnu > 0) ? nChan : 0 ) ;
2376 int izincr = ((dnu > 0) ? 1 : -1 ) ;
2377 int ichan = ((olddnu > 0) ? 0 : oldsize ) ;
2378 int iye = ((olddnu > 0) ? oldsize : 0 ) ;
2379 int iyincr = ((olddnu > 0) ? 1 : -1 ) ;
2380 //for ( int ii = izs ; ii != ize ; ii+=izincr ){
2381 int ii = izs ;
2382 while (ii != ize) {
2383 // always zl < zr
2384 double zl = zi[ii] ;
2385 double zr = zi[ii+izincr] ;
2386 // Need to access smaller index for the new spec, flag, and tsys.
2387 // Values between zi[k] and zi[k+1] should be stored in newspec[k], etc.
2388 int i = min(ii, ii+izincr) ;
2389 //for ( int jj = ichan ; jj != iye ; jj+=iyincr ) {
2390 int jj = ichan ;
2391 while (jj != iye) {
2392 // always yl < yr
2393 double yl = yi[jj] ;
2394 double yr = yi[jj+iyincr] ;
2395 // Need to access smaller index for the original spec, flag, and tsys.
2396 // Values between yi[k] and yi[k+1] are stored in oldspec[k], etc.
2397 int j = min(jj, jj+iyincr) ;
2398 if ( yr <= zl ) {
2399 jj += iyincr ;
2400 continue ;
2401 }
2402 else if ( yl <= zl ) {
2403 if ( yr < zr ) {
2404 if (!oldflag[j]) {
2405 newspec[i] += oldspec[j] * ( yr - zl ) ;
2406 if (regridTsys) newtsys[i] += oldtsys[j] * ( yr - zl ) ;
2407 wsum += ( yr - zl ) ;
2408 count[i]++ ;
2409 }
2410 newflag[i] = newflag[i] && oldflag[j] ;
2411 }
2412 else {
2413 if (!oldflag[j]) {
2414 newspec[i] += oldspec[j] * abs(dnu) ;
2415 if (regridTsys) newtsys[i] += oldtsys[j] * abs(dnu) ;
2416 wsum += abs(dnu) ;
2417 count[i]++ ;
2418 }
2419 newflag[i] = newflag[i] && oldflag[j] ;
2420 ichan = jj ;
2421 break ;
2422 }
2423 }
2424 else if ( yl < zr ) {
2425 if ( yr <= zr ) {
2426 if (!oldflag[j]) {
2427 newspec[i] += oldspec[j] * ( yr - yl ) ;
2428 if (regridTsys) newtsys[i] += oldtsys[j] * ( yr - yl ) ;
2429 wsum += ( yr - yl ) ;
2430 count[i]++ ;
2431 }
2432 newflag[i] = newflag[i] && oldflag[j] ;
2433 }
2434 else {
2435 if (!oldflag[j]) {
2436 newspec[i] += oldspec[j] * ( zr - yl ) ;
2437 if (regridTsys) newtsys[i] += oldtsys[j] * ( zr - yl ) ;
2438 wsum += ( zr - yl ) ;
2439 count[i]++ ;
2440 }
2441 newflag[i] = newflag[i] && oldflag[j] ;
2442 ichan = jj ;
2443 break ;
2444 }
2445 }
2446 else {
2447 //ichan = jj - iyincr ;
2448 break ;
2449 }
2450 jj += iyincr ;
2451 }
2452 if ( wsum != 0.0 ) {
2453 newspec[i] /= wsum ;
2454 if (regridTsys) newtsys[i] /= wsum ;
2455 }
2456 wsum = 0.0 ;
2457 ii += izincr ;
2458 }
2459
2460 // flag out channels without data
2461 // this is tentative since there is no specific definition
2462 // on bit flag...
2463 uChar noData = 1 << 7 ;
2464 for ( Int i = 0 ; i < nChan ; i++ ) {
2465 if ( count[i] == 0 )
2466 newflag[i] = noData ;
2467 }
2468
2469 specCol_.put( irow, newspec ) ;
2470 flagsCol_.put( irow, newflag ) ;
2471 if (regridTsys) tsysCol_.put( irow, newtsys );
2472
2473 return ;
2474}
2475
2476std::vector<float> Scantable::getWeather(int whichrow) const
2477{
2478 std::vector<float> out(5);
2479 //Float temperature, pressure, humidity, windspeed, windaz;
2480 weatherTable_.getEntry(out[0], out[1], out[2], out[3], out[4],
2481 mweatheridCol_(uInt(whichrow)));
2482
2483
2484 return out;
2485}
2486
2487bool Scantable::isAllChannelsFlagged(uInt whichrow)
2488{
2489 uInt rflag;
2490 flagrowCol_.get(whichrow, rflag);
2491 if (rflag > 0)
2492 return true;
2493 uChar flag;
2494 Vector<uChar> flags;
2495 flagsCol_.get(whichrow, flags);
2496 flag = flags[0];
2497 for (uInt i = 1; i < flags.size(); ++i) {
2498 flag &= flags[i];
2499 }
2500 // return ((flag >> 7) == 1);
2501 return (flag > 0);
2502}
2503
2504std::vector<std::string> Scantable::applyBaselineTable(const std::string& bltable, const bool returnfitresult, const std::string& outbltable, const bool outbltableexists, const bool overwrite)
2505{
2506 STBaselineTable btin = STBaselineTable(bltable);
2507
2508 Vector<Bool> applyCol = btin.getApply();
2509 int nRowBl = applyCol.size();
2510 if (nRowBl != nrow()) {
2511 throw(AipsError("Scantable and bltable have different number of rows."));
2512 }
2513
2514 std::vector<std::string> res;
2515 res.clear();
2516
2517 bool outBaselineTable = ((outbltable != "") && (!outbltableexists || overwrite));
2518 bool bltableidentical = (bltable == outbltable);
2519 STBaselineTable btout = STBaselineTable(*this);
2520 ROScalarColumn<Double> tcol = ROScalarColumn<Double>(table_, "TIME");
2521 Vector<Double> timeSecCol = tcol.getColumn();
2522
2523 for (int whichrow = 0; whichrow < nRowBl; ++whichrow) {
2524 if (applyCol[whichrow]) {
2525 std::vector<float> spec = getSpectrum(whichrow);
2526
2527 std::vector<bool> mask = btin.getMask(whichrow); //use mask_bltable only
2528
2529 STBaselineFunc::FuncName ftype = btin.getFunctionName(whichrow);
2530 std::vector<int> fpar = btin.getFuncParam(whichrow);
2531 std::vector<float> params;
2532 float rms;
2533 std::vector<float> resfit = doApplyBaselineTable(spec, mask, ftype, fpar, params, rms);
2534 setSpectrum(resfit, whichrow);
2535
2536 if (returnfitresult) {
2537 res.push_back(packFittingResults(whichrow, params, rms));
2538 }
2539
2540 if (outBaselineTable) {
2541 if (outbltableexists) {
2542 if (overwrite) {
2543 if (bltableidentical) {
2544 btin.setresult(uInt(whichrow), Vector<Float>(params), Float(rms));
2545 } else {
2546 btout.setresult(uInt(whichrow), Vector<Float>(params), Float(rms));
2547 }
2548 }
2549 } else {
2550 btout.appenddata(getScan(whichrow), getCycle(whichrow), getBeam(whichrow),
2551 getIF(whichrow), getPol(whichrow), 0, timeSecCol[whichrow],
2552 true, ftype, fpar, std::vector<float>(),
2553 getMaskListFromMask(mask), params, rms, spec.size(),
2554 3.0, 0, 0.0, 0, std::vector<int>());
2555 }
2556 }
2557 }
2558 }
2559
2560 if (outBaselineTable) {
2561 if (bltableidentical) {
2562 btin.save(outbltable);
2563 } else {
2564 btout.save(outbltable);
2565 }
2566 }
2567
2568 return res;
2569}
2570
2571std::vector<std::string> Scantable::subBaseline(const std::vector<std::string>& blInfoList, const bool returnfitresult, const std::string& outbltable, const bool outbltableexists, const bool overwrite)
2572{
2573 int nRowBl = blInfoList.size();
2574 int nRowSt = nrow();
2575
2576 std::vector<std::string> res;
2577 res.clear();
2578
2579 bool outBaselineTable = ((outbltable != "") && (!outbltableexists || overwrite));
2580 if ((outbltable != "") && outbltableexists && !overwrite) {
2581 throw(AipsError("Cannot overwrite bltable. Set overwrite=True."));
2582 }
2583
2584 STBaselineTable bt = STBaselineTable(*this);
2585 ROScalarColumn<Double> tcol = ROScalarColumn<Double>(table_, "TIME");
2586 Vector<Double> timeSecCol = tcol.getColumn();
2587
2588 if (outBaselineTable && !outbltableexists) {
2589 for (int i = 0; i < nRowSt; ++i) {
2590 bt.appendbasedata(getScan(i), getCycle(i), getBeam(i), getIF(i), getPol(i),
2591 0, timeSecCol[i]);
2592 bt.setApply(i, false);
2593 }
2594 }
2595
2596 for (int i = 0; i < nRowBl; ++i) {
2597 int irow;
2598 STBaselineFunc::FuncName ftype;
2599 std::vector<bool> mask;
2600 std::vector<int> fpar;
2601 float clipth;
2602 int clipn;
2603 bool uself;
2604 float lfth;
2605 std::vector<int> lfedge;
2606 int lfavg;
2607 parseBlInfo(blInfoList[i], irow, ftype, fpar, mask, clipth, clipn, uself, lfth, lfedge, lfavg);
2608
2609 if (irow < nRowSt) {
2610 std::vector<float> spec = getSpectrum(irow);
2611 std::vector<float> params;
2612 float rms;
2613 std::vector<bool> finalmask;
2614
2615 std::vector<float> resfit = doSubtractBaseline(spec, mask, ftype, fpar, params, rms, finalmask, clipth, clipn, uself, irow, lfth, lfedge, lfavg);
2616 setSpectrum(resfit, irow);
2617
2618 if (returnfitresult) {
2619 res.push_back(packFittingResults(irow, params, rms));
2620 }
2621
2622 if (outBaselineTable) {
2623 Vector<Int> fparam(fpar.size());
2624 for (uInt j = 0; j < fparam.size(); ++j) {
2625 fparam[j] = (Int)fpar[j];
2626 }
2627
2628 bt.setdata(uInt(irow),
2629 uInt(getScan(irow)), uInt(getCycle(irow)),
2630 uInt(getBeam(irow)), uInt(getIF(irow)), uInt(getPol(irow)),
2631 uInt(0), timeSecCol[irow], Bool(true), ftype, fparam,
2632 Vector<Float>(), getMaskListFromMask(finalmask), Vector<Float>(params),
2633 Float(rms), uInt(spec.size()), Float(clipth), uInt(clipn),
2634 Float(0.0), uInt(0), Vector<uInt>());
2635 }
2636
2637 }
2638 }
2639
2640 if (outBaselineTable) {
2641 bt.save(outbltable);
2642 }
2643
2644 return res;
2645}
2646
2647std::vector<float> Scantable::doApplyBaselineTable(std::vector<float>& spec,
2648 std::vector<bool>& mask,
2649 const STBaselineFunc::FuncName ftype,
2650 std::vector<int>& fpar,
2651 std::vector<float>& params,
2652 float&rms)
2653{
2654 std::vector<bool> finalmask;
2655 std::vector<int> lfedge;
2656 return doSubtractBaseline(spec, mask, ftype, fpar, params, rms, finalmask, 0.0, 0, false, 0, 0.0, lfedge, 0);
2657}
2658
2659std::vector<float> Scantable::doSubtractBaseline(std::vector<float>& spec,
2660 std::vector<bool>& mask,
2661 const STBaselineFunc::FuncName ftype,
2662 std::vector<int>& fpar,
2663 std::vector<float>& params,
2664 float&rms,
2665 std::vector<bool>& finalmask,
2666 float clipth,
2667 int clipn,
2668 bool uself,
2669 int irow,
2670 float lfth,
2671 std::vector<int>& lfedge,
2672 int lfavg)
2673{
2674 if (uself) {
2675 STLineFinder lineFinder = STLineFinder();
2676 initLineFinder(lfedge, lfth, lfavg, lineFinder);
2677 std::vector<int> currentEdge;
2678 mask = getCompositeChanMask(irow, mask, lfedge, currentEdge, lineFinder);
2679 }
2680
2681 std::vector<float> res;
2682 if (ftype == STBaselineFunc::Polynomial) {
2683 res = doPolynomialFitting(spec, mask, fpar[0], params, rms, finalmask, clipth, clipn);
2684 } else if (ftype == STBaselineFunc::Chebyshev) {
2685 res = doChebyshevFitting(spec, mask, fpar[0], params, rms, finalmask, clipth, clipn);
2686 } else if (ftype == STBaselineFunc::CSpline) {
2687 if (fpar.size() > 1) { // reading from baseline table in which pieceEdges are already calculated and stored.
2688 res = doCubicSplineFitting(spec, mask, fpar, params, rms, finalmask, clipth, clipn);
2689 } else { // usual cspline fitting by giving nPiece only. fpar will be replaced with pieceEdges.
2690 res = doCubicSplineFitting(spec, mask, fpar[0], fpar, params, rms, finalmask, clipth, clipn);
2691 }
2692 } else if (ftype == STBaselineFunc::Sinusoid) {
2693 res = doSinusoidFitting(spec, mask, fpar, params, rms, finalmask, clipth, clipn);
2694 }
2695
2696 return res;
2697}
2698
2699std::string Scantable::packFittingResults(const int irow, const std::vector<float>& params, const float rms)
2700{
2701 // returned value: "irow:params[0],params[1],..,params[n-1]:rms"
2702 ostringstream os;
2703 os << irow << ':';
2704 for (uInt i = 0; i < params.size(); ++i) {
2705 if (i > 0) {
2706 os << ',';
2707 }
2708 os << params[i];
2709 }
2710 os << ':' << rms;
2711
2712 return os.str();
2713}
2714
2715void Scantable::parseBlInfo(const std::string& blInfo, int& irow, STBaselineFunc::FuncName& ftype, std::vector<int>& fpar, std::vector<bool>& mask, float& thresClip, int& nIterClip, bool& useLineFinder, float& thresLF, std::vector<int>& edgeLF, int& avgLF)
2716{
2717 // The baseline info to be parsed must be column-delimited string like
2718 // "0:chebyshev:5:3,5,169,174,485,487" where the elements are
2719 // row number, funcType, funcOrder, maskList, clipThreshold, clipNIter,
2720 // useLineFinder, lfThreshold, lfEdge and lfChanAvgLimit.
2721
2722 std::vector<string> res = splitToStringList(blInfo, ':');
2723 if (res.size() < 4) {
2724 throw(AipsError("baseline info has bad format")) ;
2725 }
2726
2727 string ftype0, fpar0, masklist0, uself0, edge0;
2728 std::vector<int> masklist;
2729
2730 stringstream ss;
2731 ss << res[0];
2732 ss >> irow;
2733 ss.clear(); ss.str("");
2734
2735 ss << res[1];
2736 ss >> ftype0;
2737 if (ftype0 == "poly") {
2738 ftype = STBaselineFunc::Polynomial;
2739 } else if (ftype0 == "cspline") {
2740 ftype = STBaselineFunc::CSpline;
2741 } else if (ftype0 == "sinusoid") {
2742 ftype = STBaselineFunc::Sinusoid;
2743 } else if (ftype0 == "chebyshev") {
2744 ftype = STBaselineFunc::Chebyshev;
2745 } else {
2746 throw(AipsError("invalid function type."));
2747 }
2748 ss.clear(); ss.str("");
2749
2750 ss << res[2];
2751 ss >> fpar0;
2752 fpar = splitToIntList(fpar0, ',');
2753 ss.clear(); ss.str("");
2754
2755 ss << res[3];
2756 ss >> masklist0;
2757 mask = getMaskFromMaskList(nchan(getIF(irow)), splitToIntList(masklist0, ','));
2758
2759 ss << res[4];
2760 ss >> thresClip;
2761 ss.clear(); ss.str("");
2762
2763 ss << res[5];
2764 ss >> nIterClip;
2765 ss.clear(); ss.str("");
2766
2767 ss << res[6];
2768 ss >> uself0;
2769 if (uself0 == "true") {
2770 useLineFinder = true;
2771 } else {
2772 useLineFinder = false;
2773 }
2774 ss.clear(); ss.str("");
2775
2776 if (useLineFinder) {
2777 ss << res[7];
2778 ss >> thresLF;
2779 ss.clear(); ss.str("");
2780
2781 ss << res[8];
2782 ss >> edge0;
2783 edgeLF = splitToIntList(edge0, ',');
2784 ss.clear(); ss.str("");
2785
2786 ss << res[9];
2787 ss >> avgLF;
2788 ss.clear(); ss.str("");
2789 }
2790
2791}
2792
2793std::vector<int> Scantable::splitToIntList(const std::string& s, const char delim)
2794{
2795 istringstream iss(s);
2796 string tmp;
2797 int tmpi;
2798 std::vector<int> res;
2799 stringstream ss;
2800 while (getline(iss, tmp, delim)) {
2801 ss << tmp;
2802 ss >> tmpi;
2803 res.push_back(tmpi);
2804 ss.clear(); ss.str("");
2805 }
2806
2807 return res;
2808}
2809
2810std::vector<string> Scantable::splitToStringList(const std::string& s, const char delim)
2811{
2812 istringstream iss(s);
2813 std::string tmp;
2814 std::vector<string> res;
2815 while (getline(iss, tmp, delim)) {
2816 res.push_back(tmp);
2817 }
2818
2819 return res;
2820}
2821
2822std::vector<bool> Scantable::getMaskFromMaskList(const int nchan, const std::vector<int>& masklist)
2823{
2824 if (masklist.size() % 2 != 0) {
2825 throw(AipsError("masklist must have even number of elements."));
2826 }
2827
2828 std::vector<bool> res(nchan);
2829
2830 for (int i = 0; i < nchan; ++i) {
2831 res[i] = false;
2832 }
2833 for (uInt j = 0; j < masklist.size(); j += 2) {
2834 for (int i = masklist[j]; i <= masklist[j+1]; ++i) {
2835 res[i] = true;
2836 }
2837 }
2838
2839 return res;
2840}
2841
2842Vector<uInt> Scantable::getMaskListFromMask(const std::vector<bool>& mask)
2843{
2844 std::vector<int> masklist;
2845 masklist.clear();
2846
2847 for (uInt i = 0; i < mask.size(); ++i) {
2848 if (mask[i]) {
2849 if ((i == 0)||(i == mask.size()-1)) {
2850 masklist.push_back(i);
2851 } else {
2852 if ((mask[i])&&(!mask[i-1])) {
2853 masklist.push_back(i);
2854 }
2855 if ((mask[i])&&(!mask[i+1])) {
2856 masklist.push_back(i);
2857 }
2858 }
2859 }
2860 }
2861
2862 Vector<uInt> res(masklist.size());
2863 for (uInt i = 0; i < masklist.size(); ++i) {
2864 res[i] = (uInt)masklist[i];
2865 }
2866
2867 return res;
2868}
2869
2870void Scantable::initialiseBaselining(const std::string& blfile,
2871 ofstream& ofs,
2872 const bool outLogger,
2873 bool& outTextFile,
2874 bool& csvFormat,
2875 String& coordInfo,
2876 bool& hasSameNchan,
2877 const std::string& progressInfo,
2878 bool& showProgress,
2879 int& minNRow,
2880 Vector<Double>& timeSecCol)
2881{
2882 csvFormat = false;
2883 outTextFile = false;
2884
2885 if (blfile != "") {
2886 csvFormat = (blfile.substr(0, 1) == "T");
2887 ofs.open(blfile.substr(1).c_str(), ios::out | ios::app);
2888 if (ofs) outTextFile = true;
2889 }
2890
2891 coordInfo = "";
2892 hasSameNchan = true;
2893
2894 if (outLogger || outTextFile) {
2895 coordInfo = getCoordInfo()[0];
2896 if (coordInfo == "") coordInfo = "channel";
2897 hasSameNchan = hasSameNchanOverIFs();
2898 }
2899
2900 parseProgressInfo(progressInfo, showProgress, minNRow);
2901
2902 ROScalarColumn<Double> tcol = ROScalarColumn<Double>(table_, "TIME");
2903 timeSecCol = tcol.getColumn();
2904}
2905
2906void Scantable::finaliseBaselining(const bool outBaselineTable,
2907 STBaselineTable* pbt,
2908 const string& bltable,
2909 const bool outTextFile,
2910 ofstream& ofs)
2911{
2912 if (outBaselineTable) {
2913 pbt->save(bltable);
2914 }
2915
2916 if (outTextFile) ofs.close();
2917}
2918
2919void Scantable::initLineFinder(const std::vector<int>& edge,
2920 const float threshold,
2921 const int chanAvgLimit,
2922 STLineFinder& lineFinder)
2923{
2924 if ((edge.size() > 2) && (edge.size() < getIFNos().size()*2)) {
2925 throw(AipsError("Length of edge element info is less than that of IFs"));
2926 }
2927
2928 lineFinder.setOptions(threshold, 3, chanAvgLimit);
2929}
2930
2931void Scantable::polyBaseline(const std::vector<bool>& mask, int order,
2932 float thresClip, int nIterClip,
2933 bool getResidual,
2934 const std::string& progressInfo,
2935 const bool outLogger, const std::string& blfile,
2936 const std::string& bltable)
2937{
2938 /****
2939 double TimeStart = mathutil::gettimeofday_sec();
2940 ****/
2941
2942 try {
2943 ofstream ofs;
2944 String coordInfo;
2945 bool hasSameNchan, outTextFile, csvFormat, showProgress;
2946 int minNRow;
2947 int nRow = nrow();
2948 std::vector<bool> chanMask, finalChanMask;
2949 float rms;
2950 bool outBaselineTable = (bltable != "");
2951 STBaselineTable bt = STBaselineTable(*this);
2952 Vector<Double> timeSecCol;
2953
2954 initialiseBaselining(blfile, ofs, outLogger, outTextFile, csvFormat,
2955 coordInfo, hasSameNchan,
2956 progressInfo, showProgress, minNRow,
2957 timeSecCol);
2958
2959 std::vector<int> nChanNos;
2960 std::vector<std::vector<std::vector<double> > > modelReservoir;
2961 modelReservoir = getPolynomialModelReservoir(order,
2962 &Scantable::getNormalPolynomial,
2963 nChanNos);
2964
2965 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
2966 std::vector<float> sp = getSpectrum(whichrow);
2967 chanMask = getCompositeChanMask(whichrow, mask);
2968
2969 std::vector<float> params;
2970 int nClipped = 0;
2971 std::vector<float> res = doLeastSquareFitting(sp, chanMask,
2972 modelReservoir[getIdxOfNchan(sp.size(), nChanNos)],
2973 params, rms, finalChanMask,
2974 nClipped, thresClip, nIterClip, getResidual);
2975
2976 if (outBaselineTable) {
2977 bt.appenddata(getScan(whichrow), getCycle(whichrow), getBeam(whichrow),
2978 getIF(whichrow), getPol(whichrow), 0, timeSecCol[whichrow],
2979 true, STBaselineFunc::Polynomial, order, std::vector<float>(),
2980 getMaskListFromMask(finalChanMask), params, rms, sp.size(),
2981 thresClip, nIterClip, 0.0, 0, std::vector<int>());
2982 } else {
2983 setSpectrum(res, whichrow);
2984 }
2985
2986 outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow,
2987 coordInfo, hasSameNchan, ofs, "polyBaseline()",
2988 params, nClipped);
2989 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
2990 }
2991
2992 finaliseBaselining(outBaselineTable, &bt, bltable, outTextFile, ofs);
2993
2994 } catch (...) {
2995 throw;
2996 }
2997
2998 /****
2999 double TimeEnd = mathutil::gettimeofday_sec();
3000 double elapse1 = TimeEnd - TimeStart;
3001 std::cout << "poly-new : " << elapse1 << " (sec.)" << endl;
3002 ****/
3003}
3004
3005void Scantable::autoPolyBaseline(const std::vector<bool>& mask, int order,
3006 float thresClip, int nIterClip,
3007 const std::vector<int>& edge,
3008 float threshold, int chanAvgLimit,
3009 bool getResidual,
3010 const std::string& progressInfo,
3011 const bool outLogger, const std::string& blfile,
3012 const std::string& bltable)
3013{
3014 try {
3015 ofstream ofs;
3016 String coordInfo;
3017 bool hasSameNchan, outTextFile, csvFormat, showProgress;
3018 int minNRow;
3019 int nRow = nrow();
3020 std::vector<bool> chanMask, finalChanMask;
3021 float rms;
3022 bool outBaselineTable = (bltable != "");
3023 STBaselineTable bt = STBaselineTable(*this);
3024 Vector<Double> timeSecCol;
3025 STLineFinder lineFinder = STLineFinder();
3026
3027 initialiseBaselining(blfile, ofs, outLogger, outTextFile, csvFormat,
3028 coordInfo, hasSameNchan,
3029 progressInfo, showProgress, minNRow,
3030 timeSecCol);
3031
3032 initLineFinder(edge, threshold, chanAvgLimit, lineFinder);
3033
3034 std::vector<int> nChanNos;
3035 std::vector<std::vector<std::vector<double> > > modelReservoir;
3036 modelReservoir = getPolynomialModelReservoir(order,
3037 &Scantable::getNormalPolynomial,
3038 nChanNos);
3039
3040 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
3041 std::vector<float> sp = getSpectrum(whichrow);
3042 std::vector<int> currentEdge;
3043 chanMask = getCompositeChanMask(whichrow, mask, edge, currentEdge, lineFinder);
3044
3045 std::vector<float> params;
3046 int nClipped = 0;
3047 std::vector<float> res = doLeastSquareFitting(sp, chanMask,
3048 modelReservoir[getIdxOfNchan(sp.size(), nChanNos)],
3049 params, rms, finalChanMask,
3050 nClipped, thresClip, nIterClip, getResidual);
3051
3052 if (outBaselineTable) {
3053 bt.appenddata(getScan(whichrow), getCycle(whichrow), getBeam(whichrow),
3054 getIF(whichrow), getPol(whichrow), 0, timeSecCol[whichrow],
3055 true, STBaselineFunc::Polynomial, order, std::vector<float>(),
3056 getMaskListFromMask(finalChanMask), params, rms, sp.size(),
3057 thresClip, nIterClip, threshold, chanAvgLimit, currentEdge);
3058 } else {
3059 setSpectrum(res, whichrow);
3060 }
3061
3062 outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow,
3063 coordInfo, hasSameNchan, ofs, "autoPolyBaseline()",
3064 params, nClipped);
3065 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
3066 }
3067
3068 finaliseBaselining(outBaselineTable, &bt, bltable, outTextFile, ofs);
3069
3070 } catch (...) {
3071 throw;
3072 }
3073}
3074
3075void Scantable::chebyshevBaseline(const std::vector<bool>& mask, int order,
3076 float thresClip, int nIterClip,
3077 bool getResidual,
3078 const std::string& progressInfo,
3079 const bool outLogger, const std::string& blfile,
3080 const std::string& bltable)
3081{
3082 /*
3083 double TimeStart = mathutil::gettimeofday_sec();
3084 */
3085
3086 try {
3087 ofstream ofs;
3088 String coordInfo;
3089 bool hasSameNchan, outTextFile, csvFormat, showProgress;
3090 int minNRow;
3091 int nRow = nrow();
3092 std::vector<bool> chanMask, finalChanMask;
3093 float rms;
3094 bool outBaselineTable = (bltable != "");
3095 STBaselineTable bt = STBaselineTable(*this);
3096 Vector<Double> timeSecCol;
3097
3098 initialiseBaselining(blfile, ofs, outLogger, outTextFile, csvFormat,
3099 coordInfo, hasSameNchan,
3100 progressInfo, showProgress, minNRow,
3101 timeSecCol);
3102
3103 std::vector<int> nChanNos;
3104 std::vector<std::vector<std::vector<double> > > modelReservoir;
3105 modelReservoir = getPolynomialModelReservoir(order,
3106 &Scantable::getChebyshevPolynomial,
3107 nChanNos);
3108
3109 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
3110 std::vector<float> sp = getSpectrum(whichrow);
3111 chanMask = getCompositeChanMask(whichrow, mask);
3112
3113 std::vector<float> params;
3114 int nClipped = 0;
3115 std::vector<float> res = doLeastSquareFitting(sp, chanMask,
3116 modelReservoir[getIdxOfNchan(sp.size(), nChanNos)],
3117 params, rms, finalChanMask,
3118 nClipped, thresClip, nIterClip, getResidual);
3119
3120 if (outBaselineTable) {
3121 bt.appenddata(getScan(whichrow), getCycle(whichrow), getBeam(whichrow),
3122 getIF(whichrow), getPol(whichrow), 0, timeSecCol[whichrow],
3123 true, STBaselineFunc::Chebyshev, order, std::vector<float>(),
3124 getMaskListFromMask(finalChanMask), params, rms, sp.size(),
3125 thresClip, nIterClip, 0.0, 0, std::vector<int>());
3126 } else {
3127 setSpectrum(res, whichrow);
3128 }
3129
3130 outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow,
3131 coordInfo, hasSameNchan, ofs, "chebyshevBaseline()",
3132 params, nClipped);
3133 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
3134 }
3135
3136 finaliseBaselining(outBaselineTable, &bt, bltable, outTextFile, ofs);
3137
3138 } catch (...) {
3139 throw;
3140 }
3141
3142 /*
3143 double TimeEnd = mathutil::gettimeofday_sec();
3144 double elapse1 = TimeEnd - TimeStart;
3145 std::cout << "cheby : " << elapse1 << " (sec.)" << endl;
3146 */
3147}
3148
3149void Scantable::autoChebyshevBaseline(const std::vector<bool>& mask, int order,
3150 float thresClip, int nIterClip,
3151 const std::vector<int>& edge,
3152 float threshold, int chanAvgLimit,
3153 bool getResidual,
3154 const std::string& progressInfo,
3155 const bool outLogger, const std::string& blfile,
3156 const std::string& bltable)
3157{
3158 try {
3159 ofstream ofs;
3160 String coordInfo;
3161 bool hasSameNchan, outTextFile, csvFormat, showProgress;
3162 int minNRow;
3163 int nRow = nrow();
3164 std::vector<bool> chanMask, finalChanMask;
3165 float rms;
3166 bool outBaselineTable = (bltable != "");
3167 STBaselineTable bt = STBaselineTable(*this);
3168 Vector<Double> timeSecCol;
3169 STLineFinder lineFinder = STLineFinder();
3170
3171 initialiseBaselining(blfile, ofs, outLogger, outTextFile, csvFormat,
3172 coordInfo, hasSameNchan,
3173 progressInfo, showProgress, minNRow,
3174 timeSecCol);
3175
3176 initLineFinder(edge, threshold, chanAvgLimit, lineFinder);
3177
3178 std::vector<int> nChanNos;
3179 std::vector<std::vector<std::vector<double> > > modelReservoir;
3180 modelReservoir = getPolynomialModelReservoir(order,
3181 &Scantable::getChebyshevPolynomial,
3182 nChanNos);
3183
3184 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
3185 std::vector<float> sp = getSpectrum(whichrow);
3186 std::vector<int> currentEdge;
3187 chanMask = getCompositeChanMask(whichrow, mask, edge, currentEdge, lineFinder);
3188
3189 std::vector<float> params;
3190 int nClipped = 0;
3191 std::vector<float> res = doLeastSquareFitting(sp, chanMask,
3192 modelReservoir[getIdxOfNchan(sp.size(), nChanNos)],
3193 params, rms, finalChanMask,
3194 nClipped, thresClip, nIterClip, getResidual);
3195
3196 if (outBaselineTable) {
3197 bt.appenddata(getScan(whichrow), getCycle(whichrow), getBeam(whichrow),
3198 getIF(whichrow), getPol(whichrow), 0, timeSecCol[whichrow],
3199 true, STBaselineFunc::Chebyshev, order, std::vector<float>(),
3200 getMaskListFromMask(finalChanMask), params, rms, sp.size(),
3201 thresClip, nIterClip, threshold, chanAvgLimit, currentEdge);
3202 } else {
3203 setSpectrum(res, whichrow);
3204 }
3205
3206 outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow,
3207 coordInfo, hasSameNchan, ofs, "autoChebyshevBaseline()",
3208 params, nClipped);
3209 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
3210 }
3211
3212 finaliseBaselining(outBaselineTable, &bt, bltable, outTextFile, ofs);
3213
3214 } catch (...) {
3215 throw;
3216 }
3217}
3218
3219double Scantable::calculateModelSelectionCriteria(const std::string& valname,
3220 const std::string& blfunc,
3221 int order,
3222 const std::vector<bool>& inMask,
3223 int whichrow,
3224 bool useLineFinder,
3225 const std::vector<int>& edge,
3226 float threshold,
3227 int chanAvgLimit)
3228{
3229 std::vector<float> sp = getSpectrum(whichrow);
3230 std::vector<bool> chanMask;
3231 chanMask.clear();
3232
3233 if (useLineFinder) {
3234 STLineFinder lineFinder = STLineFinder();
3235 initLineFinder(edge, threshold, chanAvgLimit, lineFinder);
3236 std::vector<int> currentEdge;
3237 chanMask = getCompositeChanMask(whichrow, inMask, edge, currentEdge, lineFinder);
3238 } else {
3239 chanMask = getCompositeChanMask(whichrow, inMask);
3240 }
3241
3242 return doCalculateModelSelectionCriteria(valname, sp, chanMask, blfunc, order);
3243}
3244
3245double Scantable::doCalculateModelSelectionCriteria(const std::string& valname, const std::vector<float>& spec, const std::vector<bool>& mask, const std::string& blfunc, int order)
3246{
3247 int nparam;
3248 std::vector<float> params;
3249 std::vector<bool> finalChanMask;
3250 float rms;
3251 int nClipped = 0;
3252 std::vector<float> res;
3253 if (blfunc == "poly") {
3254 nparam = order + 1;
3255 res = doPolynomialFitting(spec, mask, order, params, rms, finalChanMask, nClipped);
3256 } else if (blfunc == "chebyshev") {
3257 nparam = order + 1;
3258 res = doChebyshevFitting(spec, mask, order, params, rms, finalChanMask, nClipped);
3259 } else if (blfunc == "cspline") {
3260 std::vector<int> pieceEdges;//(order+1); //order = npiece
3261 nparam = order + 3;
3262 res = doCubicSplineFitting(spec, mask, order, false, pieceEdges, params, rms, finalChanMask, nClipped);
3263 } else if (blfunc == "sinusoid") {
3264 std::vector<int> nWaves;
3265 nWaves.clear();
3266 for (int i = 0; i <= order; ++i) {
3267 nWaves.push_back(i);
3268 }
3269 nparam = 2*order + 1; // order = nwave
3270 res = doSinusoidFitting(spec, mask, nWaves, params, rms, finalChanMask, nClipped);
3271 } else {
3272 throw(AipsError("blfunc must be poly, chebyshev, cspline or sinusoid."));
3273 }
3274
3275 double msq = 0.0;
3276 int nusedchan = 0;
3277 int nChan = res.size();
3278 for (int i = 0; i < nChan; ++i) {
3279 if (mask[i]) {
3280 msq += (double)res[i]*(double)res[i];
3281 nusedchan++;
3282 }
3283 }
3284 if (nusedchan == 0) {
3285 throw(AipsError("all channels masked."));
3286 }
3287 msq /= (double)nusedchan;
3288
3289 nparam++; //add 1 for sigma of Gaussian distribution
3290 const double PI = 6.0 * asin(0.5); // PI (= 3.141592653...)
3291
3292 if (valname.find("aic") == 0) {
3293 // Original Akaike Information Criterion (AIC)
3294 double aic = nusedchan * (log(2.0 * PI * msq) + 1.0) + 2.0 * nparam;
3295
3296 // Corrected AIC by Sugiura(1978) (AICc)
3297 if (valname == "aicc") {
3298 if (nusedchan - nparam - 1 <= 0) {
3299 throw(AipsError("channel size is too small to calculate AICc."));
3300 }
3301 aic += 2.0*nparam*(nparam + 1)/(double)(nusedchan - nparam - 1);
3302 }
3303
3304 return aic;
3305
3306 } else if (valname == "bic") {
3307 // Bayesian Information Criterion (BIC)
3308 double bic = nusedchan * log(msq) + nparam * log((double)nusedchan);
3309 return bic;
3310
3311 } else if (valname == "gcv") {
3312 // Generalised Cross Validation
3313 double x = 1.0 - (double)nparam / (double)nusedchan;
3314 double gcv = msq / (x * x);
3315 return gcv;
3316
3317 } else {
3318 throw(AipsError("valname must be aic, aicc, bic or gcv."));
3319 }
3320}
3321
3322double Scantable::getNormalPolynomial(int n, double x) {
3323 if (n == 0) {
3324 return 1.0;
3325 } else if (n > 0) {
3326 double res = 1.0;
3327 for (int i = 0; i < n; ++i) {
3328 res *= x;
3329 }
3330 return res;
3331 } else {
3332 if (x == 0.0) {
3333 throw(AipsError("infinity result: x=0 given for negative power."));
3334 } else {
3335 return pow(x, (double)n);
3336 }
3337 }
3338}
3339
3340double Scantable::getChebyshevPolynomial(int n, double x) {
3341 if ((x < -1.0)||(x > 1.0)) {
3342 throw(AipsError("out of definition range (-1 <= x <= 1)."));
3343 } else if (x == 1.0) {
3344 return 1.0;
3345 } else if (x == 0.0) {
3346 double res;
3347 if (n%2 == 0) {
3348 if (n%4 == 0) {
3349 res = 1.0;
3350 } else {
3351 res = -1.0;
3352 }
3353 } else {
3354 res = 0.0;
3355 }
3356 return res;
3357 } else if (x == -1.0) {
3358 double res = (n%2 == 0 ? 1.0 : -1.0);
3359 return res;
3360 } else if (n < 0) {
3361 throw(AipsError("the order must be zero or positive."));
3362 } else if (n == 0) {
3363 return 1.0;
3364 } else if (n == 1) {
3365 return x;
3366 } else {
3367 double res = 0.0;
3368 for (int m = 0; m <= n/2; ++m) {
3369 double c = 1.0;
3370 if (m > 0) {
3371 for (int i = 1; i <= m; ++i) {
3372 c *= (double)(n-2*m+i)/(double)i;
3373 }
3374 }
3375 res += (m%2 == 0 ? 1.0 : -1.0)*(double)n/(double)(n-m)*pow(2.0*x, (double)(n-2*m))/2.0*c;
3376 }
3377 return res;
3378 }
3379}
3380
3381std::vector<float> Scantable::doPolynomialFitting(const std::vector<float>& data,
3382 const std::vector<bool>& mask,
3383 int order,
3384 std::vector<float>& params,
3385 float& rms,
3386 std::vector<bool>& finalmask,
3387 float clipth,
3388 int clipn)
3389{
3390 int nClipped = 0;
3391 return doPolynomialFitting(data, mask, order, params, rms, finalmask, nClipped, clipth, clipn);
3392}
3393
3394std::vector<float> Scantable::doPolynomialFitting(const std::vector<float>& data,
3395 const std::vector<bool>& mask,
3396 int order,
3397 std::vector<float>& params,
3398 float& rms,
3399 std::vector<bool>& finalMask,
3400 int& nClipped,
3401 float thresClip,
3402 int nIterClip,
3403 bool getResidual)
3404{
3405 return doLeastSquareFitting(data, mask,
3406 getPolynomialModel(order, data.size(), &Scantable::getNormalPolynomial),
3407 params, rms, finalMask,
3408 nClipped, thresClip, nIterClip,
3409 getResidual);
3410}
3411
3412std::vector<float> Scantable::doChebyshevFitting(const std::vector<float>& data,
3413 const std::vector<bool>& mask,
3414 int order,
3415 std::vector<float>& params,
3416 float& rms,
3417 std::vector<bool>& finalmask,
3418 float clipth,
3419 int clipn)
3420{
3421 int nClipped = 0;
3422 return doChebyshevFitting(data, mask, order, params, rms, finalmask, nClipped, clipth, clipn);
3423}
3424
3425std::vector<float> Scantable::doChebyshevFitting(const std::vector<float>& data,
3426 const std::vector<bool>& mask,
3427 int order,
3428 std::vector<float>& params,
3429 float& rms,
3430 std::vector<bool>& finalMask,
3431 int& nClipped,
3432 float thresClip,
3433 int nIterClip,
3434 bool getResidual)
3435{
3436 return doLeastSquareFitting(data, mask,
3437 getPolynomialModel(order, data.size(), &Scantable::getChebyshevPolynomial),
3438 params, rms, finalMask,
3439 nClipped, thresClip, nIterClip,
3440 getResidual);
3441}
3442
3443std::vector<std::vector<double> > Scantable::getPolynomialModel(int order, int nchan, double (Scantable::*pfunc)(int, double))
3444{
3445 // model : contains model values for computing the least-square matrix.
3446 // model.size() is nmodel and model[*].size() is nchan.
3447 // Each model element are as follows:
3448 //
3449 // (for normal polynomials)
3450 // model[0] = {1.0, 1.0, 1.0, ..., 1.0},
3451 // model[1] = {0.0, 1.0, 2.0, ..., (nchan-1)}
3452 // model[n-1] = ...,
3453 // model[n] = {0.0^n, 1.0^n, 2.0^n, ..., (nchan-1)^n}
3454 // where (0 <= n <= order)
3455 //
3456 // (for Chebyshev polynomials)
3457 // model[0] = {T0(-1), T0(2/(nchan-1)-1), T0(4/(nchan-1)-1), ..., T0(1)},
3458 // model[n-1] = ...,
3459 // model[n] = {Tn(-1), Tn(2/(nchan-1)-1), Tn(4/(nchan-1)-1), ..., Tn(1)}
3460 // where (0 <= n <= order),
3461
3462 int nmodel = order + 1;
3463 std::vector<std::vector<double> > model(nmodel, std::vector<double>(nchan));
3464
3465 double stretch, shift;
3466 if (pfunc == &Scantable::getChebyshevPolynomial) {
3467 stretch = 2.0/(double)(nchan - 1);
3468 shift = -1.0;
3469 } else {
3470 stretch = 1.0;
3471 shift = 0.0;
3472 }
3473
3474 for (int i = 0; i < nmodel; ++i) {
3475 for (int j = 0; j < nchan; ++j) {
3476 model[i][j] = (this->*pfunc)(i, stretch*(double)j + shift);
3477 }
3478 }
3479
3480 return model;
3481}
3482
3483std::vector<std::vector<std::vector<double> > > Scantable::getPolynomialModelReservoir(int order,
3484 double (Scantable::*pfunc)(int, double),
3485 std::vector<int>& nChanNos)
3486{
3487 std::vector<std::vector<std::vector<double> > > res;
3488 res.clear();
3489 nChanNos.clear();
3490
3491 std::vector<uint> ifNos = getIFNos();
3492 for (uint i = 0; i < ifNos.size(); ++i) {
3493 int currNchan = nchan(ifNos[i]);
3494 bool hasDifferentNchan = (i == 0);
3495 for (uint j = 0; j < i; ++j) {
3496 if (currNchan != nchan(ifNos[j])) {
3497 hasDifferentNchan = true;
3498 break;
3499 }
3500 }
3501 if (hasDifferentNchan) {
3502 res.push_back(getPolynomialModel(order, currNchan, pfunc));
3503 nChanNos.push_back(currNchan);
3504 }
3505 }
3506
3507 return res;
3508}
3509
3510std::vector<float> Scantable::doLeastSquareFitting(const std::vector<float>& data,
3511 const std::vector<bool>& mask,
3512 const std::vector<std::vector<double> >& model,
3513 std::vector<float>& params,
3514 float& rms,
3515 std::vector<bool>& finalMask,
3516 int& nClipped,
3517 float thresClip,
3518 int nIterClip,
3519 bool getResidual)
3520{
3521 int nDOF = model.size();
3522 int nChan = data.size();
3523
3524 if (nDOF == 0) {
3525 throw(AipsError("no model data given"));
3526 }
3527 if (nChan < 2) {
3528 throw(AipsError("data size is too few"));
3529 }
3530 if (nChan != (int)mask.size()) {
3531 throw(AipsError("data and mask sizes are not identical"));
3532 }
3533 for (int i = 0; i < nDOF; ++i) {
3534 if (nChan != (int)model[i].size()) {
3535 throw(AipsError("data and model sizes are not identical"));
3536 }
3537 }
3538
3539 params.clear();
3540 params.resize(nDOF);
3541
3542 finalMask.clear();
3543 finalMask.resize(nChan);
3544
3545 std::vector<int> maskArray(nChan);
3546 int j = 0;
3547 for (int i = 0; i < nChan; ++i) {
3548 maskArray[i] = mask[i] ? 1 : 0;
3549 if (mask[i]) {
3550 j++;
3551 }
3552 finalMask[i] = mask[i];
3553 }
3554
3555 int initNData = j;
3556 int nData = initNData;
3557
3558 std::vector<double> z1(nChan), r1(nChan), residual(nChan);
3559 for (int i = 0; i < nChan; ++i) {
3560 z1[i] = (double)data[i];
3561 r1[i] = 0.0;
3562 residual[i] = 0.0;
3563 }
3564
3565 for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
3566 // xMatrix : horizontal concatenation of
3567 // the least-sq. matrix (left) and an
3568 // identity matrix (right).
3569 // the right part is used to calculate the inverse matrix of the left part.
3570 double xMatrix[nDOF][2*nDOF];
3571 double zMatrix[nDOF];
3572 for (int i = 0; i < nDOF; ++i) {
3573 for (int j = 0; j < 2*nDOF; ++j) {
3574 xMatrix[i][j] = 0.0;
3575 }
3576 xMatrix[i][nDOF+i] = 1.0;
3577 zMatrix[i] = 0.0;
3578 }
3579
3580 int nUseData = 0;
3581 for (int k = 0; k < nChan; ++k) {
3582 if (maskArray[k] == 0) continue;
3583
3584 for (int i = 0; i < nDOF; ++i) {
3585 for (int j = i; j < nDOF; ++j) {
3586 xMatrix[i][j] += model[i][k] * model[j][k];
3587 }
3588 zMatrix[i] += z1[k] * model[i][k];
3589 }
3590
3591 nUseData++;
3592 }
3593
3594 if (nUseData < 1) {
3595 throw(AipsError("all channels clipped or masked. can't execute fitting anymore."));
3596 }
3597
3598 for (int i = 0; i < nDOF; ++i) {
3599 for (int j = 0; j < i; ++j) {
3600 xMatrix[i][j] = xMatrix[j][i];
3601 }
3602 }
3603
3604 std::vector<double> invDiag(nDOF);
3605 for (int i = 0; i < nDOF; ++i) {
3606 invDiag[i] = 1.0 / xMatrix[i][i];
3607 for (int j = 0; j < nDOF; ++j) {
3608 xMatrix[i][j] *= invDiag[i];
3609 }
3610 }
3611
3612 for (int k = 0; k < nDOF; ++k) {
3613 for (int i = 0; i < nDOF; ++i) {
3614 if (i != k) {
3615 double factor1 = xMatrix[k][k];
3616 double invfactor1 = 1.0 / factor1;
3617 double factor2 = xMatrix[i][k];
3618 for (int j = k; j < 2*nDOF; ++j) {
3619 xMatrix[i][j] *= factor1;
3620 xMatrix[i][j] -= xMatrix[k][j]*factor2;
3621 xMatrix[i][j] *= invfactor1;
3622 }
3623 }
3624 }
3625 double invXDiag = 1.0 / xMatrix[k][k];
3626 for (int j = k; j < 2*nDOF; ++j) {
3627 xMatrix[k][j] *= invXDiag;
3628 }
3629 }
3630
3631 for (int i = 0; i < nDOF; ++i) {
3632 for (int j = 0; j < nDOF; ++j) {
3633 xMatrix[i][nDOF+j] *= invDiag[j];
3634 }
3635 }
3636 //compute a vector y in which coefficients of the best-fit
3637 //model functions are stored.
3638 //in case of polynomials, y consists of (a0,a1,a2,...)
3639 //where ai is the coefficient of the term x^i.
3640 //in case of sinusoids, y consists of (a0,s1,c1,s2,c2,...)
3641 //where a0 is constant term and s* and c* are of sine
3642 //and cosine functions, respectively.
3643 std::vector<double> y(nDOF);
3644 for (int i = 0; i < nDOF; ++i) {
3645 y[i] = 0.0;
3646 for (int j = 0; j < nDOF; ++j) {
3647 y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
3648 }
3649 params[i] = (float)y[i];
3650 }
3651
3652 for (int i = 0; i < nChan; ++i) {
3653 r1[i] = y[0];
3654 for (int j = 1; j < nDOF; ++j) {
3655 r1[i] += y[j]*model[j][i];
3656 }
3657 residual[i] = z1[i] - r1[i];
3658 }
3659
3660 double stdDev = 0.0;
3661 for (int i = 0; i < nChan; ++i) {
3662 stdDev += residual[i]*residual[i]*(double)maskArray[i];
3663 }
3664 stdDev = sqrt(stdDev/(double)nData);
3665 rms = (float)stdDev;
3666
3667 if ((nClip == nIterClip) || (thresClip <= 0.0)) {
3668 break;
3669 } else {
3670
3671 double thres = stdDev * thresClip;
3672 int newNData = 0;
3673 for (int i = 0; i < nChan; ++i) {
3674 if (abs(residual[i]) >= thres) {
3675 maskArray[i] = 0;
3676 finalMask[i] = false;
3677 }
3678 if (maskArray[i] > 0) {
3679 newNData++;
3680 }
3681 }
3682 if (newNData == nData) {
3683 break; //no more flag to add. iteration stops.
3684 } else {
3685 nData = newNData;
3686 }
3687
3688 }
3689 }
3690
3691 nClipped = initNData - nData;
3692
3693 std::vector<float> result(nChan);
3694 if (getResidual) {
3695 for (int i = 0; i < nChan; ++i) {
3696 result[i] = (float)residual[i];
3697 }
3698 } else {
3699 for (int i = 0; i < nChan; ++i) {
3700 result[i] = (float)r1[i];
3701 }
3702 }
3703
3704 return result;
3705}
3706
3707void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece,
3708 float thresClip, int nIterClip,
3709 bool getResidual,
3710 const std::string& progressInfo,
3711 const bool outLogger, const std::string& blfile,
3712 const std::string& bltable)
3713{
3714 /****
3715 double TimeStart = mathutil::gettimeofday_sec();
3716 ****/
3717
3718 try {
3719 ofstream ofs;
3720 String coordInfo;
3721 bool hasSameNchan, outTextFile, csvFormat, showProgress;
3722 int minNRow;
3723 int nRow = nrow();
3724 std::vector<bool> chanMask, finalChanMask;
3725 float rms;
3726 bool outBaselineTable = (bltable != "");
3727 STBaselineTable bt = STBaselineTable(*this);
3728 Vector<Double> timeSecCol;
3729
3730 initialiseBaselining(blfile, ofs, outLogger, outTextFile, csvFormat,
3731 coordInfo, hasSameNchan,
3732 progressInfo, showProgress, minNRow,
3733 timeSecCol);
3734
3735 std::vector<int> nChanNos;
3736 std::vector<std::vector<std::vector<double> > > modelReservoir;
3737 modelReservoir = getPolynomialModelReservoir(3,
3738 &Scantable::getNormalPolynomial,
3739 nChanNos);
3740
3741 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
3742 std::vector<float> sp = getSpectrum(whichrow);
3743 chanMask = getCompositeChanMask(whichrow, mask);
3744
3745 std::vector<int> pieceEdges;
3746 std::vector<float> params;
3747 int nClipped = 0;
3748 std::vector<float> res = doCubicSplineLeastSquareFitting(sp, chanMask,
3749 modelReservoir[getIdxOfNchan(sp.size(), nChanNos)],
3750 nPiece, false, pieceEdges, params, rms, finalChanMask,
3751 nClipped, thresClip, nIterClip, getResidual);
3752
3753 if (outBaselineTable) {
3754 bt.appenddata(getScan(whichrow), getCycle(whichrow), getBeam(whichrow),
3755 getIF(whichrow), getPol(whichrow), 0, timeSecCol[whichrow],
3756 true, STBaselineFunc::CSpline, pieceEdges, std::vector<float>(),
3757 getMaskListFromMask(finalChanMask), params, rms, sp.size(),
3758 thresClip, nIterClip, 0.0, 0, std::vector<int>());
3759 } else {
3760 setSpectrum(res, whichrow);
3761 }
3762
3763 outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow,
3764 coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()",
3765 pieceEdges, params, nClipped);
3766 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
3767 }
3768
3769 finaliseBaselining(outBaselineTable, &bt, bltable, outTextFile, ofs);
3770
3771 } catch (...) {
3772 throw;
3773 }
3774
3775 /****
3776 double TimeEnd = mathutil::gettimeofday_sec();
3777 double elapse1 = TimeEnd - TimeStart;
3778 std::cout << "cspline-new : " << elapse1 << " (sec.)" << endl;
3779 ****/
3780}
3781
3782void Scantable::autoCubicSplineBaseline(const std::vector<bool>& mask, int nPiece,
3783 float thresClip, int nIterClip,
3784 const std::vector<int>& edge,
3785 float threshold, int chanAvgLimit,
3786 bool getResidual,
3787 const std::string& progressInfo,
3788 const bool outLogger, const std::string& blfile,
3789 const std::string& bltable)
3790{
3791 try {
3792 ofstream ofs;
3793 String coordInfo;
3794 bool hasSameNchan, outTextFile, csvFormat, showProgress;
3795 int minNRow;
3796 int nRow = nrow();
3797 std::vector<bool> chanMask, finalChanMask;
3798 float rms;
3799 bool outBaselineTable = (bltable != "");
3800 STBaselineTable bt = STBaselineTable(*this);
3801 Vector<Double> timeSecCol;
3802 STLineFinder lineFinder = STLineFinder();
3803
3804 initialiseBaselining(blfile, ofs, outLogger, outTextFile, csvFormat,
3805 coordInfo, hasSameNchan,
3806 progressInfo, showProgress, minNRow,
3807 timeSecCol);
3808
3809 initLineFinder(edge, threshold, chanAvgLimit, lineFinder);
3810
3811 std::vector<int> nChanNos;
3812 std::vector<std::vector<std::vector<double> > > modelReservoir;
3813 modelReservoir = getPolynomialModelReservoir(3,
3814 &Scantable::getNormalPolynomial,
3815 nChanNos);
3816
3817 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
3818 std::vector<float> sp = getSpectrum(whichrow);
3819 std::vector<int> currentEdge;
3820 chanMask = getCompositeChanMask(whichrow, mask, edge, currentEdge, lineFinder);
3821
3822 std::vector<int> pieceEdges;
3823 std::vector<float> params;
3824 int nClipped = 0;
3825 std::vector<float> res = doCubicSplineLeastSquareFitting(sp, chanMask,
3826 modelReservoir[getIdxOfNchan(sp.size(), nChanNos)],
3827 nPiece, false, pieceEdges, params, rms, finalChanMask,
3828 nClipped, thresClip, nIterClip, getResidual);
3829
3830 if (outBaselineTable) {
3831 bt.appenddata(getScan(whichrow), getCycle(whichrow), getBeam(whichrow),
3832 getIF(whichrow), getPol(whichrow), 0, timeSecCol[whichrow],
3833 true, STBaselineFunc::CSpline, pieceEdges, std::vector<float>(),
3834 getMaskListFromMask(finalChanMask), params, rms, sp.size(),
3835 thresClip, nIterClip, threshold, chanAvgLimit, currentEdge);
3836 } else {
3837 setSpectrum(res, whichrow);
3838 }
3839
3840 outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow,
3841 coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()",
3842 pieceEdges, params, nClipped);
3843 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
3844 }
3845
3846 finaliseBaselining(outBaselineTable, &bt, bltable, outTextFile, ofs);
3847
3848 } catch (...) {
3849 throw;
3850 }
3851}
3852
3853std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data,
3854 const std::vector<bool>& mask,
3855 std::vector<int>& idxEdge,
3856 std::vector<float>& params,
3857 float& rms,
3858 std::vector<bool>& finalmask,
3859 float clipth,
3860 int clipn)
3861{
3862 int nClipped = 0;
3863 return doCubicSplineFitting(data, mask, idxEdge.size()-1, true, idxEdge, params, rms, finalmask, nClipped, clipth, clipn);
3864}
3865
3866std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data,
3867 const std::vector<bool>& mask,
3868 int nPiece,
3869 std::vector<int>& idxEdge,
3870 std::vector<float>& params,
3871 float& rms,
3872 std::vector<bool>& finalmask,
3873 float clipth,
3874 int clipn)
3875{
3876 int nClipped = 0;
3877 return doCubicSplineFitting(data, mask, nPiece, false, idxEdge, params, rms, finalmask, nClipped, clipth, clipn);
3878}
3879
3880std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data,
3881 const std::vector<bool>& mask,
3882 int nPiece,
3883 bool useGivenPieceBoundary,
3884 std::vector<int>& idxEdge,
3885 std::vector<float>& params,
3886 float& rms,
3887 std::vector<bool>& finalMask,
3888 int& nClipped,
3889 float thresClip,
3890 int nIterClip,
3891 bool getResidual)
3892{
3893 return doCubicSplineLeastSquareFitting(data, mask,
3894 getPolynomialModel(3, data.size(), &Scantable::getNormalPolynomial),
3895 nPiece, useGivenPieceBoundary, idxEdge,
3896 params, rms, finalMask,
3897 nClipped, thresClip, nIterClip,
3898 getResidual);
3899}
3900
3901std::vector<float> Scantable::doCubicSplineLeastSquareFitting(const std::vector<float>& data,
3902 const std::vector<bool>& mask,
3903 const std::vector<std::vector<double> >& model,
3904 int nPiece,
3905 bool useGivenPieceBoundary,
3906 std::vector<int>& idxEdge,
3907 std::vector<float>& params,
3908 float& rms,
3909 std::vector<bool>& finalMask,
3910 int& nClipped,
3911 float thresClip,
3912 int nIterClip,
3913 bool getResidual)
3914{
3915 int nDOF = nPiece + 3; //number of independent parameters to solve, namely, 4+(nPiece-1).
3916 int nModel = model.size();
3917 int nChan = data.size();
3918
3919 if (nModel != 4) {
3920 throw(AipsError("model size must be 4."));
3921 }
3922 if (nPiece < 1) {
3923 throw(AipsError("number of the sections must be one or more"));
3924 }
3925 if (nChan < 2*nPiece) {
3926 throw(AipsError("data size is too few"));
3927 }
3928 if (nChan != (int)mask.size()) {
3929 throw(AipsError("data and mask sizes are not identical"));
3930 }
3931 for (int i = 0; i < nModel; ++i) {
3932 if (nChan != (int)model[i].size()) {
3933 throw(AipsError("data and model sizes are not identical"));
3934 }
3935 }
3936
3937 params.clear();
3938 params.resize(nPiece*nModel);
3939
3940 finalMask.clear();
3941 finalMask.resize(nChan);
3942
3943 std::vector<int> maskArray(nChan);
3944 std::vector<int> x(nChan);
3945 int j = 0;
3946 for (int i = 0; i < nChan; ++i) {
3947 maskArray[i] = mask[i] ? 1 : 0;
3948 if (mask[i]) {
3949 x[j] = i;
3950 j++;
3951 }
3952 finalMask[i] = mask[i];
3953 }
3954
3955 int initNData = j;
3956 int nData = initNData;
3957
3958 if (initNData < nPiece) {
3959 throw(AipsError("too few non-flagged channels"));
3960 }
3961
3962 int nElement = (int)(floor(floor((double)(initNData/nPiece))+0.5));
3963 std::vector<double> invEdge(nPiece-1);
3964
3965 if (useGivenPieceBoundary) {
3966 if ((int)idxEdge.size() != nPiece+1) {
3967 throw(AipsError("pieceEdge.size() must be equal to nPiece+1."));
3968 }
3969 } else {
3970 idxEdge.clear();
3971 idxEdge.resize(nPiece+1);
3972 idxEdge[0] = x[0];
3973 }
3974 for (int i = 1; i < nPiece; ++i) {
3975 int valX = x[nElement*i];
3976 if (!useGivenPieceBoundary) {
3977 idxEdge[i] = valX;
3978 }
3979 invEdge[i-1] = 1.0/(double)valX;
3980 }
3981 if (!useGivenPieceBoundary) {
3982 idxEdge[nPiece] = x[initNData-1]+1;
3983 }
3984
3985 std::vector<double> z1(nChan), r1(nChan), residual(nChan);
3986 for (int i = 0; i < nChan; ++i) {
3987 z1[i] = (double)data[i];
3988 r1[i] = 0.0;
3989 residual[i] = 0.0;
3990 }
3991
3992 for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
3993 // xMatrix : horizontal concatenation of
3994 // the least-sq. matrix (left) and an
3995 // identity matrix (right).
3996 // the right part is used to calculate the inverse matrix of the left part.
3997
3998 double xMatrix[nDOF][2*nDOF];
3999 double zMatrix[nDOF];
4000 for (int i = 0; i < nDOF; ++i) {
4001 for (int j = 0; j < 2*nDOF; ++j) {
4002 xMatrix[i][j] = 0.0;
4003 }
4004 xMatrix[i][nDOF+i] = 1.0;
4005 zMatrix[i] = 0.0;
4006 }
4007
4008 for (int n = 0; n < nPiece; ++n) {
4009 int nUseDataInPiece = 0;
4010 for (int k = idxEdge[n]; k < idxEdge[n+1]; ++k) {
4011
4012 if (maskArray[k] == 0) continue;
4013
4014 for (int i = 0; i < nModel; ++i) {
4015 for (int j = i; j < nModel; ++j) {
4016 xMatrix[i][j] += model[i][k] * model[j][k];
4017 }
4018 zMatrix[i] += z1[k] * model[i][k];
4019 }
4020
4021 for (int i = 0; i < n; ++i) {
4022 double q = 1.0 - model[1][k]*invEdge[i];
4023 q = q*q*q;
4024 for (int j = 0; j < nModel; ++j) {
4025 xMatrix[j][i+nModel] += q * model[j][k];
4026 }
4027 for (int j = 0; j < i; ++j) {
4028 double r = 1.0 - model[1][k]*invEdge[j];
4029 r = r*r*r;
4030 xMatrix[j+nModel][i+nModel] += r*q;
4031 }
4032 xMatrix[i+nModel][i+nModel] += q*q;
4033 zMatrix[i+nModel] += q*z1[k];
4034 }
4035
4036 nUseDataInPiece++;
4037 }
4038
4039 if (nUseDataInPiece < 1) {
4040 std::vector<string> suffixOfPieceNumber(4);
4041 suffixOfPieceNumber[0] = "th";
4042 suffixOfPieceNumber[1] = "st";
4043 suffixOfPieceNumber[2] = "nd";
4044 suffixOfPieceNumber[3] = "rd";
4045 int idxNoDataPiece = (n % 10 <= 3) ? n : 0;
4046 ostringstream oss;
4047 oss << "all channels clipped or masked in " << n << suffixOfPieceNumber[idxNoDataPiece];
4048 oss << " piece of the spectrum. can't execute fitting anymore.";
4049 throw(AipsError(String(oss)));
4050 }
4051 }
4052
4053 for (int i = 0; i < nDOF; ++i) {
4054 for (int j = 0; j < i; ++j) {
4055 xMatrix[i][j] = xMatrix[j][i];
4056 }
4057 }
4058
4059 std::vector<double> invDiag(nDOF);
4060 for (int i = 0; i < nDOF; ++i) {
4061 invDiag[i] = 1.0 / xMatrix[i][i];
4062 for (int j = 0; j < nDOF; ++j) {
4063 xMatrix[i][j] *= invDiag[i];
4064 }
4065 }
4066
4067 for (int k = 0; k < nDOF; ++k) {
4068 for (int i = 0; i < nDOF; ++i) {
4069 if (i != k) {
4070 double factor1 = xMatrix[k][k];
4071 double invfactor1 = 1.0 / factor1;
4072 double factor2 = xMatrix[i][k];
4073 for (int j = k; j < 2*nDOF; ++j) {
4074 xMatrix[i][j] *= factor1;
4075 xMatrix[i][j] -= xMatrix[k][j]*factor2;
4076 xMatrix[i][j] *= invfactor1;
4077 }
4078 }
4079 }
4080 double invXDiag = 1.0 / xMatrix[k][k];
4081 for (int j = k; j < 2*nDOF; ++j) {
4082 xMatrix[k][j] *= invXDiag;
4083 }
4084 }
4085
4086 for (int i = 0; i < nDOF; ++i) {
4087 for (int j = 0; j < nDOF; ++j) {
4088 xMatrix[i][nDOF+j] *= invDiag[j];
4089 }
4090 }
4091
4092 //compute a vector y which consists of the coefficients of the best-fit spline curves
4093 //(a0,a1,a2,a3(,b3,c3,...)), namely, the ones for the leftmost piece and the ones of
4094 //cubic terms for the other pieces (in case nPiece>1).
4095 std::vector<double> y(nDOF);
4096 for (int i = 0; i < nDOF; ++i) {
4097 y[i] = 0.0;
4098 for (int j = 0; j < nDOF; ++j) {
4099 y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
4100 }
4101 }
4102
4103 std::vector<double> a(nModel);
4104 for (int i = 0; i < nModel; ++i) {
4105 a[i] = y[i];
4106 }
4107
4108 int j = 0;
4109 for (int n = 0; n < nPiece; ++n) {
4110 for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
4111 r1[i] = 0.0;
4112 for (int j = 0; j < nModel; ++j) {
4113 r1[i] += a[j] * model[j][i];
4114 }
4115 }
4116 for (int i = 0; i < nModel; ++i) {
4117 params[j+i] = a[i];
4118 }
4119 j += nModel;
4120
4121 if (n == nPiece-1) break;
4122
4123 double d = y[n+nModel];
4124 double iE = invEdge[n];
4125 a[0] += d;
4126 a[1] -= 3.0 * d * iE;
4127 a[2] += 3.0 * d * iE * iE;
4128 a[3] -= d * iE * iE * iE;
4129 }
4130
4131 //subtract constant value for masked regions at the edge of spectrum
4132 if (idxEdge[0] > 0) {
4133 int n = idxEdge[0];
4134 for (int i = 0; i < idxEdge[0]; ++i) {
4135 //--cubic extrapolate--
4136 //r1[i] = params[0] + params[1]*x1[i] + params[2]*x2[i] + params[3]*x3[i];
4137 //--linear extrapolate--
4138 //r1[i] = (r1[n+1] - r1[n])/(x1[n+1] - x1[n])*(x1[i] - x1[n]) + r1[n];
4139 //--constant--
4140 r1[i] = r1[n];
4141 }
4142 }
4143
4144 if (idxEdge[nPiece] < nChan) {
4145 int n = idxEdge[nPiece]-1;
4146 for (int i = idxEdge[nPiece]; i < nChan; ++i) {
4147 //--cubic extrapolate--
4148 //int m = 4*(nPiece-1);
4149 //r1[i] = params[m] + params[m+1]*x1[i] + params[m+2]*x2[i] + params[m+3]*x3[i];
4150 //--linear extrapolate--
4151 //r1[i] = (r1[n-1] - r1[n])/(x1[n-1] - x1[n])*(x1[i] - x1[n]) + r1[n];
4152 //--constant--
4153 r1[i] = r1[n];
4154 }
4155 }
4156
4157 for (int i = 0; i < nChan; ++i) {
4158 residual[i] = z1[i] - r1[i];
4159 }
4160
4161 double stdDev = 0.0;
4162 for (int i = 0; i < nChan; ++i) {
4163 stdDev += residual[i]*residual[i]*(double)maskArray[i];
4164 }
4165 stdDev = sqrt(stdDev/(double)nData);
4166 rms = (float)stdDev;
4167
4168 if ((nClip == nIterClip) || (thresClip <= 0.0)) {
4169 break;
4170 } else {
4171
4172 double thres = stdDev * thresClip;
4173 int newNData = 0;
4174 for (int i = 0; i < nChan; ++i) {
4175 if (abs(residual[i]) >= thres) {
4176 maskArray[i] = 0;
4177 finalMask[i] = false;
4178 }
4179 if (maskArray[i] > 0) {
4180 newNData++;
4181 }
4182 }
4183 if (newNData == nData) {
4184 break; //no more flag to add. iteration stops.
4185 } else {
4186 nData = newNData;
4187 }
4188
4189 }
4190 }
4191
4192 nClipped = initNData - nData;
4193
4194 std::vector<float> result(nChan);
4195 if (getResidual) {
4196 for (int i = 0; i < nChan; ++i) {
4197 result[i] = (float)residual[i];
4198 }
4199 } else {
4200 for (int i = 0; i < nChan; ++i) {
4201 result[i] = (float)r1[i];
4202 }
4203 }
4204
4205 return result;
4206}
4207
4208std::vector<int> Scantable::selectWaveNumbers(const std::vector<int>& addNWaves,
4209 const std::vector<int>& rejectNWaves)
4210{
4211 std::vector<bool> chanMask;
4212 std::string fftMethod;
4213 std::string fftThresh;
4214
4215 return selectWaveNumbers(0, chanMask, false, fftMethod, fftThresh, addNWaves, rejectNWaves);
4216}
4217
4218std::vector<int> Scantable::selectWaveNumbers(const int whichrow,
4219 const std::vector<bool>& chanMask,
4220 const bool applyFFT,
4221 const std::string& fftMethod,
4222 const std::string& fftThresh,
4223 const std::vector<int>& addNWaves,
4224 const std::vector<int>& rejectNWaves)
4225{
4226 std::vector<int> nWaves;
4227 nWaves.clear();
4228
4229 if (applyFFT) {
4230 string fftThAttr;
4231 float fftThSigma;
4232 int fftThTop;
4233 parseFFTThresholdInfo(fftThresh, fftThAttr, fftThSigma, fftThTop);
4234 doSelectWaveNumbers(whichrow, chanMask, fftMethod, fftThSigma, fftThTop, fftThAttr, nWaves);
4235 }
4236
4237 addAuxWaveNumbers(whichrow, addNWaves, rejectNWaves, nWaves);
4238
4239 return nWaves;
4240}
4241
4242int Scantable::getIdxOfNchan(const int nChan, const std::vector<int>& nChanNos)
4243{
4244 int idx = -1;
4245 for (uint i = 0; i < nChanNos.size(); ++i) {
4246 if (nChan == nChanNos[i]) {
4247 idx = i;
4248 break;
4249 }
4250 }
4251
4252 if (idx < 0) {
4253 throw(AipsError("nChan not found in nChhanNos."));
4254 }
4255
4256 return idx;
4257}
4258
4259void Scantable::parseFFTInfo(const std::string& fftInfo, bool& applyFFT, std::string& fftMethod, std::string& fftThresh)
4260{
4261 istringstream iss(fftInfo);
4262 std::string tmp;
4263 std::vector<string> res;
4264 while (getline(iss, tmp, ',')) {
4265 res.push_back(tmp);
4266 }
4267 if (res.size() < 3) {
4268 throw(AipsError("wrong value in 'fftinfo' parameter")) ;
4269 }
4270 applyFFT = (res[0] == "true");
4271 fftMethod = res[1];
4272 fftThresh = res[2];
4273}
4274
4275void Scantable::parseFFTThresholdInfo(const std::string& fftThresh, std::string& fftThAttr, float& fftThSigma, int& fftThTop)
4276{
4277 uInt idxSigma = fftThresh.find("sigma");
4278 uInt idxTop = fftThresh.find("top");
4279
4280 if (idxSigma == fftThresh.size() - 5) {
4281 std::istringstream is(fftThresh.substr(0, fftThresh.size() - 5));
4282 is >> fftThSigma;
4283 fftThAttr = "sigma";
4284 } else if (idxTop == 0) {
4285 std::istringstream is(fftThresh.substr(3));
4286 is >> fftThTop;
4287 fftThAttr = "top";
4288 } else {
4289 bool isNumber = true;
4290 for (uInt i = 0; i < fftThresh.size()-1; ++i) {
4291 char ch = (fftThresh.substr(i, 1).c_str())[0];
4292 if (!(isdigit(ch) || (fftThresh.substr(i, 1) == "."))) {
4293 isNumber = false;
4294 break;
4295 }
4296 }
4297 if (isNumber) {
4298 std::istringstream is(fftThresh);
4299 is >> fftThSigma;
4300 fftThAttr = "sigma";
4301 } else {
4302 throw(AipsError("fftthresh has a wrong value"));
4303 }
4304 }
4305}
4306
4307void Scantable::doSelectWaveNumbers(const int whichrow, const std::vector<bool>& chanMask, const std::string& fftMethod, const float fftThSigma, const int fftThTop, const std::string& fftThAttr, std::vector<int>& nWaves)
4308{
4309 std::vector<float> fspec;
4310 if (fftMethod == "fft") {
4311 fspec = execFFT(whichrow, chanMask, false, true);
4312 //} else if (fftMethod == "lsp") {
4313 // fspec = lombScarglePeriodogram(whichrow);
4314 }
4315
4316 if (fftThAttr == "sigma") {
4317 float mean = 0.0;
4318 float mean2 = 0.0;
4319 for (uInt i = 0; i < fspec.size(); ++i) {
4320 mean += fspec[i];
4321 mean2 += fspec[i]*fspec[i];
4322 }
4323 mean /= float(fspec.size());
4324 mean2 /= float(fspec.size());
4325 float thres = mean + fftThSigma * float(sqrt(mean2 - mean*mean));
4326
4327 for (uInt i = 0; i < fspec.size(); ++i) {
4328 if (fspec[i] >= thres) {
4329 nWaves.push_back(i);
4330 }
4331 }
4332
4333 } else if (fftThAttr == "top") {
4334 for (int i = 0; i < fftThTop; ++i) {
4335 float max = 0.0;
4336 int maxIdx = 0;
4337 for (uInt j = 0; j < fspec.size(); ++j) {
4338 if (fspec[j] > max) {
4339 max = fspec[j];
4340 maxIdx = j;
4341 }
4342 }
4343 nWaves.push_back(maxIdx);
4344 fspec[maxIdx] = 0.0;
4345 }
4346
4347 }
4348
4349 if (nWaves.size() > 1) {
4350 sort(nWaves.begin(), nWaves.end());
4351 }
4352}
4353
4354void Scantable::addAuxWaveNumbers(const int whichrow, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
4355{
4356 std::vector<int> tempAddNWaves, tempRejectNWaves;
4357 tempAddNWaves.clear();
4358 tempRejectNWaves.clear();
4359
4360 for (uInt i = 0; i < addNWaves.size(); ++i) {
4361 tempAddNWaves.push_back(addNWaves[i]);
4362 }
4363 if ((tempAddNWaves.size() == 2) && (tempAddNWaves[1] == -999)) {
4364 setWaveNumberListUptoNyquistFreq(whichrow, tempAddNWaves);
4365 }
4366
4367 for (uInt i = 0; i < rejectNWaves.size(); ++i) {
4368 tempRejectNWaves.push_back(rejectNWaves[i]);
4369 }
4370 if ((tempRejectNWaves.size() == 2) && (tempRejectNWaves[1] == -999)) {
4371 setWaveNumberListUptoNyquistFreq(whichrow, tempRejectNWaves);
4372 }
4373
4374 for (uInt i = 0; i < tempAddNWaves.size(); ++i) {
4375 bool found = false;
4376 for (uInt j = 0; j < nWaves.size(); ++j) {
4377 if (nWaves[j] == tempAddNWaves[i]) {
4378 found = true;
4379 break;
4380 }
4381 }
4382 if (!found) nWaves.push_back(tempAddNWaves[i]);
4383 }
4384
4385 for (uInt i = 0; i < tempRejectNWaves.size(); ++i) {
4386 for (std::vector<int>::iterator j = nWaves.begin(); j != nWaves.end(); ) {
4387 if (*j == tempRejectNWaves[i]) {
4388 j = nWaves.erase(j);
4389 } else {
4390 ++j;
4391 }
4392 }
4393 }
4394
4395 if (nWaves.size() > 1) {
4396 sort(nWaves.begin(), nWaves.end());
4397 unique(nWaves.begin(), nWaves.end());
4398 }
4399}
4400
4401void Scantable::setWaveNumberListUptoNyquistFreq(const int whichrow, std::vector<int>& nWaves)
4402{
4403 int val = nWaves[0];
4404 int nyquistFreq = nchan(getIF(whichrow))/2+1;
4405 nWaves.clear();
4406 if (val > nyquistFreq) { // for safety, at least nWaves contains a constant; CAS-3759
4407 nWaves.push_back(0);
4408 }
4409 while (val <= nyquistFreq) {
4410 nWaves.push_back(val);
4411 val++;
4412 }
4413}
4414
4415void Scantable::sinusoidBaseline(const std::vector<bool>& mask, const std::string& fftInfo,
4416 const std::vector<int>& addNWaves,
4417 const std::vector<int>& rejectNWaves,
4418 float thresClip, int nIterClip,
4419 bool getResidual,
4420 const std::string& progressInfo,
4421 const bool outLogger, const std::string& blfile,
4422 const std::string& bltable)
4423{
4424 /****
4425 double TimeStart = mathutil::gettimeofday_sec();
4426 ****/
4427
4428 try {
4429 ofstream ofs;
4430 String coordInfo;
4431 bool hasSameNchan, outTextFile, csvFormat, showProgress;
4432 int minNRow;
4433 int nRow = nrow();
4434 std::vector<bool> chanMask, finalChanMask;
4435 float rms;
4436 bool outBaselineTable = (bltable != "");
4437 STBaselineTable bt = STBaselineTable(*this);
4438 Vector<Double> timeSecCol;
4439
4440 initialiseBaselining(blfile, ofs, outLogger, outTextFile, csvFormat,
4441 coordInfo, hasSameNchan,
4442 progressInfo, showProgress, minNRow,
4443 timeSecCol);
4444
4445 bool applyFFT;
4446 std::string fftMethod, fftThresh;
4447 parseFFTInfo(fftInfo, applyFFT, fftMethod, fftThresh);
4448
4449 std::vector<int> nWaves;
4450 std::vector<int> nChanNos;
4451 std::vector<std::vector<std::vector<double> > > modelReservoir;
4452 if (!applyFFT) {
4453 nWaves = selectWaveNumbers(addNWaves, rejectNWaves);
4454 modelReservoir = getSinusoidModelReservoir(nWaves, nChanNos);
4455 }
4456
4457 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
4458 std::vector<float> sp = getSpectrum(whichrow);
4459 chanMask = getCompositeChanMask(whichrow, mask);
4460 std::vector<std::vector<double> > model;
4461 if (applyFFT) {
4462 nWaves = selectWaveNumbers(whichrow, chanMask, true, fftMethod, fftThresh,
4463 addNWaves, rejectNWaves);
4464 model = getSinusoidModel(nWaves, sp.size());
4465 } else {
4466 model = modelReservoir[getIdxOfNchan(sp.size(), nChanNos)];
4467 }
4468
4469 std::vector<float> params;
4470 int nClipped = 0;
4471 std::vector<float> res = doLeastSquareFitting(sp, chanMask, model,
4472 params, rms, finalChanMask,
4473 nClipped, thresClip, nIterClip, getResidual);
4474
4475 if (outBaselineTable) {
4476 bt.appenddata(getScan(whichrow), getCycle(whichrow), getBeam(whichrow),
4477 getIF(whichrow), getPol(whichrow), 0, timeSecCol[whichrow],
4478 true, STBaselineFunc::Sinusoid, nWaves, std::vector<float>(),
4479 getMaskListFromMask(finalChanMask), params, rms, sp.size(),
4480 thresClip, nIterClip, 0.0, 0, std::vector<int>());
4481 } else {
4482 setSpectrum(res, whichrow);
4483 }
4484
4485 outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow,
4486 coordInfo, hasSameNchan, ofs, "sinusoidBaseline()",
4487 params, nClipped);
4488 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
4489 }
4490
4491 finaliseBaselining(outBaselineTable, &bt, bltable, outTextFile, ofs);
4492
4493 } catch (...) {
4494 throw;
4495 }
4496
4497 /****
4498 double TimeEnd = mathutil::gettimeofday_sec();
4499 double elapse1 = TimeEnd - TimeStart;
4500 std::cout << "sinusoid-old : " << elapse1 << " (sec.)" << endl;
4501 ****/
4502}
4503
4504void Scantable::autoSinusoidBaseline(const std::vector<bool>& mask, const std::string& fftInfo,
4505 const std::vector<int>& addNWaves,
4506 const std::vector<int>& rejectNWaves,
4507 float thresClip, int nIterClip,
4508 const std::vector<int>& edge,
4509 float threshold, int chanAvgLimit,
4510 bool getResidual,
4511 const std::string& progressInfo,
4512 const bool outLogger, const std::string& blfile,
4513 const std::string& bltable)
4514{
4515 try {
4516 ofstream ofs;
4517 String coordInfo;
4518 bool hasSameNchan, outTextFile, csvFormat, showProgress;
4519 int minNRow;
4520 int nRow = nrow();
4521 std::vector<bool> chanMask, finalChanMask;
4522 float rms;
4523 bool outBaselineTable = (bltable != "");
4524 STBaselineTable bt = STBaselineTable(*this);
4525 Vector<Double> timeSecCol;
4526 STLineFinder lineFinder = STLineFinder();
4527
4528 initialiseBaselining(blfile, ofs, outLogger, outTextFile, csvFormat,
4529 coordInfo, hasSameNchan,
4530 progressInfo, showProgress, minNRow,
4531 timeSecCol);
4532
4533 initLineFinder(edge, threshold, chanAvgLimit, lineFinder);
4534
4535 bool applyFFT;
4536 string fftMethod, fftThresh;
4537 parseFFTInfo(fftInfo, applyFFT, fftMethod, fftThresh);
4538
4539 std::vector<int> nWaves;
4540 std::vector<int> nChanNos;
4541 std::vector<std::vector<std::vector<double> > > modelReservoir;
4542 if (!applyFFT) {
4543 nWaves = selectWaveNumbers(addNWaves, rejectNWaves);
4544 modelReservoir = getSinusoidModelReservoir(nWaves, nChanNos);
4545 }
4546
4547 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
4548 std::vector<float> sp = getSpectrum(whichrow);
4549 std::vector<int> currentEdge;
4550 chanMask = getCompositeChanMask(whichrow, mask, edge, currentEdge, lineFinder);
4551 std::vector<std::vector<double> > model;
4552 if (applyFFT) {
4553 nWaves = selectWaveNumbers(whichrow, chanMask, true, fftMethod, fftThresh,
4554 addNWaves, rejectNWaves);
4555 model = getSinusoidModel(nWaves, sp.size());
4556 } else {
4557 model = modelReservoir[getIdxOfNchan(sp.size(), nChanNos)];
4558 }
4559
4560 std::vector<float> params;
4561 int nClipped = 0;
4562 std::vector<float> res = doLeastSquareFitting(sp, chanMask, model,
4563 params, rms, finalChanMask,
4564 nClipped, thresClip, nIterClip, getResidual);
4565
4566 if (outBaselineTable) {
4567 bt.appenddata(getScan(whichrow), getCycle(whichrow), getBeam(whichrow),
4568 getIF(whichrow), getPol(whichrow), 0, timeSecCol[whichrow],
4569 true, STBaselineFunc::Sinusoid, nWaves, std::vector<float>(),
4570 getMaskListFromMask(finalChanMask), params, rms, sp.size(),
4571 thresClip, nIterClip, threshold, chanAvgLimit, currentEdge);
4572 } else {
4573 setSpectrum(res, whichrow);
4574 }
4575
4576 outputFittingResult(outLogger, outTextFile, csvFormat, chanMask, whichrow,
4577 coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()",
4578 params, nClipped);
4579 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
4580 }
4581
4582 finaliseBaselining(outBaselineTable, &bt, bltable, outTextFile, ofs);
4583
4584 } catch (...) {
4585 throw;
4586 }
4587}
4588
4589std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data,
4590 const std::vector<bool>& mask,
4591 const std::vector<int>& waveNumbers,
4592 std::vector<float>& params,
4593 float& rms,
4594 std::vector<bool>& finalmask,
4595 float clipth,
4596 int clipn)
4597{
4598 int nClipped = 0;
4599 return doSinusoidFitting(data, mask, waveNumbers, params, rms, finalmask, nClipped, clipth, clipn);
4600}
4601
4602std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data,
4603 const std::vector<bool>& mask,
4604 const std::vector<int>& waveNumbers,
4605 std::vector<float>& params,
4606 float& rms,
4607 std::vector<bool>& finalMask,
4608 int& nClipped,
4609 float thresClip,
4610 int nIterClip,
4611 bool getResidual)
4612{
4613 return doLeastSquareFitting(data, mask,
4614 getSinusoidModel(waveNumbers, data.size()),
4615 params, rms, finalMask,
4616 nClipped, thresClip, nIterClip,
4617 getResidual);
4618}
4619
4620std::vector<std::vector<std::vector<double> > > Scantable::getSinusoidModelReservoir(const std::vector<int>& waveNumbers,
4621 std::vector<int>& nChanNos)
4622{
4623 std::vector<std::vector<std::vector<double> > > res;
4624 res.clear();
4625 nChanNos.clear();
4626
4627 std::vector<uint> ifNos = getIFNos();
4628 for (uint i = 0; i < ifNos.size(); ++i) {
4629 int currNchan = nchan(ifNos[i]);
4630 bool hasDifferentNchan = (i == 0);
4631 for (uint j = 0; j < i; ++j) {
4632 if (currNchan != nchan(ifNos[j])) {
4633 hasDifferentNchan = true;
4634 break;
4635 }
4636 }
4637 if (hasDifferentNchan) {
4638 res.push_back(getSinusoidModel(waveNumbers, currNchan));
4639 nChanNos.push_back(currNchan);
4640 }
4641 }
4642
4643 return res;
4644}
4645
4646std::vector<std::vector<double> > Scantable::getSinusoidModel(const std::vector<int>& waveNumbers, int nchan)
4647{
4648 // model : contains elemental values for computing the least-square matrix.
4649 // model.size() is nmodel and model[*].size() is nchan.
4650 // Each model element are as follows:
4651 // model[0] = {1.0, 1.0, 1.0, ..., 1.0},
4652 // model[2n-1] = {sin(nPI/L*x[0]), sin(nPI/L*x[1]), ..., sin(nPI/L*x[nchan])},
4653 // model[2n] = {cos(nPI/L*x[0]), cos(nPI/L*x[1]), ..., cos(nPI/L*x[nchan])},
4654 // where (1 <= n <= nMaxWavesInSW),
4655 // or,
4656 // model[2n-1] = {sin(wn[n]PI/L*x[0]), sin(wn[n]PI/L*x[1]), ..., sin(wn[n]PI/L*x[nchan])},
4657 // model[2n] = {cos(wn[n]PI/L*x[0]), cos(wn[n]PI/L*x[1]), ..., cos(wn[n]PI/L*x[nchan])},
4658 // where wn[n] denotes waveNumbers[n] (1 <= n <= waveNumbers.size()).
4659
4660 std::vector<int> nWaves; // sorted and uniqued array of wave numbers
4661 nWaves.reserve(waveNumbers.size());
4662 copy(waveNumbers.begin(), waveNumbers.end(), back_inserter(nWaves));
4663 sort(nWaves.begin(), nWaves.end());
4664 std::vector<int>::iterator end_it = unique(nWaves.begin(), nWaves.end());
4665 nWaves.erase(end_it, nWaves.end());
4666
4667 int minNWaves = nWaves[0];
4668 if (minNWaves < 0) {
4669 throw(AipsError("wave number must be positive or zero (i.e. constant)"));
4670 }
4671 bool hasConstantTerm = (minNWaves == 0);
4672 int nmodel = nWaves.size() * 2 - (hasConstantTerm ? 1 : 0); //number of parameters to solve.
4673
4674 std::vector<std::vector<double> > model(nmodel, std::vector<double>(nchan));
4675
4676 if (hasConstantTerm) {
4677 for (int j = 0; j < nchan; ++j) {
4678 model[0][j] = 1.0;
4679 }
4680 }
4681
4682 const double PI = 6.0 * asin(0.5); // PI (= 3.141592653...)
4683 double stretch0 = 2.0*PI/(double)(nchan-1);
4684
4685 for (uInt i = (hasConstantTerm ? 1 : 0); i < nWaves.size(); ++i) {
4686 int sidx = hasConstantTerm ? 2*i-1 : 2*i;
4687 int cidx = sidx + 1;
4688 double stretch = stretch0*(double)nWaves[i];
4689
4690 for (int j = 0; j < nchan; ++j) {
4691 model[sidx][j] = sin(stretch*(double)j);
4692 model[cidx][j] = cos(stretch*(double)j);
4693 }
4694 }
4695
4696 return model;
4697}
4698
4699std::vector<bool> Scantable::getCompositeChanMask(int whichrow,
4700 const std::vector<bool>& inMask)
4701{
4702 std::vector<bool> mask = getMask(whichrow);
4703 uInt maskSize = mask.size();
4704 if (inMask.size() != 0) {
4705 if (maskSize != inMask.size()) {
4706 throw(AipsError("mask sizes are not the same."));
4707 }
4708 for (uInt i = 0; i < maskSize; ++i) {
4709 mask[i] = mask[i] && inMask[i];
4710 }
4711 }
4712
4713 return mask;
4714}
4715
4716std::vector<bool> Scantable::getCompositeChanMask(int whichrow,
4717 const std::vector<bool>& inMask,
4718 const std::vector<int>& edge,
4719 std::vector<int>& currEdge,
4720 STLineFinder& lineFinder)
4721{
4722 std::vector<uint> ifNos = getIFNos();
4723 if ((edge.size() > 2) && (edge.size() < ifNos.size()*2)) {
4724 throw(AipsError("Length of edge element info is less than that of IFs"));
4725 }
4726
4727 uint idx = 0;
4728 if (edge.size() > 2) {
4729 int ifVal = getIF(whichrow);
4730 bool foundIF = false;
4731 for (uint i = 0; i < ifNos.size(); ++i) {
4732 if (ifVal == (int)ifNos[i]) {
4733 idx = 2*i;
4734 foundIF = true;
4735 break;
4736 }
4737 }
4738 if (!foundIF) {
4739 throw(AipsError("bad IF number"));
4740 }
4741 }
4742
4743 currEdge.clear();
4744 currEdge.resize(2);
4745 currEdge[0] = edge[idx];
4746 currEdge[1] = edge[idx+1];
4747
4748 lineFinder.setData(getSpectrum(whichrow));
4749 lineFinder.findLines(getCompositeChanMask(whichrow, inMask), currEdge, whichrow);
4750
4751 return lineFinder.getMask();
4752}
4753
4754/* for cspline. will be merged once cspline is available in fitter (2011/3/10 WK) */
4755void Scantable::outputFittingResult(bool outLogger,
4756 bool outTextFile,
4757 bool csvFormat,
4758 const std::vector<bool>& chanMask,
4759 int whichrow,
4760 const casa::String& coordInfo,
4761 bool hasSameNchan,
4762 ofstream& ofs,
4763 const casa::String& funcName,
4764 const std::vector<int>& edge,
4765 const std::vector<float>& params,
4766 const int nClipped)
4767{
4768 if (outLogger || outTextFile) {
4769 float rms = getRms(chanMask, whichrow);
4770 String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
4771 std::vector<bool> fixed;
4772 fixed.clear();
4773
4774 if (outLogger) {
4775 LogIO ols(LogOrigin("Scantable", funcName, WHERE));
4776 ols << formatPiecewiseBaselineParams(edge, params, fixed, rms, nClipped,
4777 masklist, whichrow, false, csvFormat) << LogIO::POST ;
4778 }
4779 if (outTextFile) {
4780 ofs << formatPiecewiseBaselineParams(edge, params, fixed, rms, nClipped,
4781 masklist, whichrow, true, csvFormat) << flush;
4782 }
4783 }
4784}
4785
4786/* for poly/chebyshev/sinusoid. */
4787void Scantable::outputFittingResult(bool outLogger,
4788 bool outTextFile,
4789 bool csvFormat,
4790 const std::vector<bool>& chanMask,
4791 int whichrow,
4792 const casa::String& coordInfo,
4793 bool hasSameNchan,
4794 ofstream& ofs,
4795 const casa::String& funcName,
4796 const std::vector<float>& params,
4797 const int nClipped)
4798{
4799 if (outLogger || outTextFile) {
4800 float rms = getRms(chanMask, whichrow);
4801 String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
4802 std::vector<bool> fixed;
4803 fixed.clear();
4804
4805 if (outLogger) {
4806 LogIO ols(LogOrigin("Scantable", funcName, WHERE));
4807 ols << formatBaselineParams(params, fixed, rms, nClipped,
4808 masklist, whichrow, false, csvFormat) << LogIO::POST ;
4809 }
4810 if (outTextFile) {
4811 ofs << formatBaselineParams(params, fixed, rms, nClipped,
4812 masklist, whichrow, true, csvFormat) << flush;
4813 }
4814 }
4815}
4816
4817void Scantable::parseProgressInfo(const std::string& progressInfo, bool& showProgress, int& minNRow)
4818{
4819 int idxDelimiter = progressInfo.find(",");
4820 if (idxDelimiter < 0) {
4821 throw(AipsError("wrong value in 'showprogress' parameter")) ;
4822 }
4823 showProgress = (progressInfo.substr(0, idxDelimiter) == "true");
4824 std::istringstream is(progressInfo.substr(idxDelimiter+1));
4825 is >> minNRow;
4826}
4827
4828void Scantable::showProgressOnTerminal(const int nProcessed, const int nTotal, const bool showProgress, const int nTotalThreshold)
4829{
4830 if (showProgress && (nTotal >= nTotalThreshold)) {
4831 int nInterval = int(floor(double(nTotal)/100.0));
4832 if (nInterval == 0) nInterval++;
4833
4834 if (nProcessed % nInterval == 0) {
4835 printf("\r"); //go to the head of line
4836 printf("\x1b[31m\x1b[1m"); //set red color, highlighted
4837 printf("[%3d%%]", (int)(100.0*(double(nProcessed+1))/(double(nTotal))) );
4838 printf("\x1b[39m\x1b[0m"); //set default attributes
4839 fflush(NULL);
4840 }
4841
4842 if (nProcessed == nTotal - 1) {
4843 printf("\r\x1b[K"); //clear
4844 fflush(NULL);
4845 }
4846
4847 }
4848}
4849
4850std::vector<float> Scantable::execFFT(const int whichrow, const std::vector<bool>& inMask, bool getRealImag, bool getAmplitudeOnly)
4851{
4852 std::vector<bool> mask = getMask(whichrow);
4853
4854 if (inMask.size() > 0) {
4855 uInt maskSize = mask.size();
4856 if (maskSize != inMask.size()) {
4857 throw(AipsError("mask sizes are not the same."));
4858 }
4859 for (uInt i = 0; i < maskSize; ++i) {
4860 mask[i] = mask[i] && inMask[i];
4861 }
4862 }
4863
4864 Vector<Float> spec = getSpectrum(whichrow);
4865 mathutil::doZeroOrderInterpolation(spec, mask);
4866
4867 FFTServer<Float,Complex> ffts;
4868 Vector<Complex> fftres;
4869 ffts.fft0(fftres, spec);
4870
4871 std::vector<float> res;
4872 float norm = float(2.0/double(spec.size()));
4873
4874 if (getRealImag) {
4875 for (uInt i = 0; i < fftres.size(); ++i) {
4876 res.push_back(real(fftres[i])*norm);
4877 res.push_back(imag(fftres[i])*norm);
4878 }
4879 } else {
4880 for (uInt i = 0; i < fftres.size(); ++i) {
4881 res.push_back(abs(fftres[i])*norm);
4882 if (!getAmplitudeOnly) res.push_back(arg(fftres[i]));
4883 }
4884 }
4885
4886 return res;
4887}
4888
4889
4890float Scantable::getRms(const std::vector<bool>& mask, int whichrow)
4891{
4892 /****
4893 double ms1TimeStart, ms1TimeEnd;
4894 double elapse1 = 0.0;
4895 ms1TimeStart = mathutil::gettimeofday_sec();
4896 ****/
4897
4898 Vector<Float> spec;
4899 specCol_.get(whichrow, spec);
4900
4901 /****
4902 ms1TimeEnd = mathutil::gettimeofday_sec();
4903 elapse1 = ms1TimeEnd - ms1TimeStart;
4904 std::cout << "rm1 : " << elapse1 << " (sec.)" << endl;
4905 ****/
4906
4907 return (float)doGetRms(mask, spec);
4908}
4909
4910double Scantable::doGetRms(const std::vector<bool>& mask, const Vector<Float>& spec)
4911{
4912 double mean = 0.0;
4913 double smean = 0.0;
4914 int n = 0;
4915 for (uInt i = 0; i < spec.nelements(); ++i) {
4916 if (mask[i]) {
4917 double val = (double)spec[i];
4918 mean += val;
4919 smean += val*val;
4920 n++;
4921 }
4922 }
4923
4924 mean /= (double)n;
4925 smean /= (double)n;
4926
4927 return sqrt(smean - mean*mean);
4928}
4929
4930std::string Scantable::formatBaselineParamsHeader(int whichrow, const std::string& masklist, bool verbose, bool csvformat) const
4931{
4932 if (verbose) {
4933 ostringstream oss;
4934
4935 if (csvformat) {
4936 oss << getScan(whichrow) << ",";
4937 oss << getBeam(whichrow) << ",";
4938 oss << getIF(whichrow) << ",";
4939 oss << getPol(whichrow) << ",";
4940 oss << getCycle(whichrow) << ",";
4941 String commaReplacedMasklist = masklist;
4942 string::size_type pos = 0;
4943 while (pos = commaReplacedMasklist.find(","), pos != string::npos) {
4944 commaReplacedMasklist.replace(pos, 1, ";");
4945 pos++;
4946 }
4947 oss << commaReplacedMasklist << ",";
4948 } else {
4949 oss << " Scan[" << getScan(whichrow) << "]";
4950 oss << " Beam[" << getBeam(whichrow) << "]";
4951 oss << " IF[" << getIF(whichrow) << "]";
4952 oss << " Pol[" << getPol(whichrow) << "]";
4953 oss << " Cycle[" << getCycle(whichrow) << "]: " << endl;
4954 oss << "Fitter range = " << masklist << endl;
4955 oss << "Baseline parameters" << endl;
4956 }
4957 oss << flush;
4958
4959 return String(oss);
4960 }
4961
4962 return "";
4963}
4964
4965std::string Scantable::formatBaselineParamsFooter(float rms, int nClipped, bool verbose, bool csvformat) const
4966{
4967 if (verbose) {
4968 ostringstream oss;
4969
4970 if (csvformat) {
4971 oss << rms << ",";
4972 if (nClipped >= 0) {
4973 oss << nClipped;
4974 }
4975 } else {
4976 oss << "Results of baseline fit" << endl;
4977 oss << " rms = " << setprecision(6) << rms << endl;
4978 if (nClipped >= 0) {
4979 oss << " Number of clipped channels = " << nClipped << endl;
4980 }
4981 for (int i = 0; i < 60; ++i) {
4982 oss << "-";
4983 }
4984 }
4985 oss << endl;
4986 oss << flush;
4987
4988 return String(oss);
4989 }
4990
4991 return "";
4992}
4993
4994std::string Scantable::formatBaselineParams(const std::vector<float>& params,
4995 const std::vector<bool>& fixed,
4996 float rms,
4997 int nClipped,
4998 const std::string& masklist,
4999 int whichrow,
5000 bool verbose,
5001 bool csvformat,
5002 int start, int count,
5003 bool resetparamid) const
5004{
5005 int nParam = (int)(params.size());
5006
5007 if (nParam < 1) {
5008 return(" Not fitted");
5009 } else {
5010
5011 ostringstream oss;
5012 oss << formatBaselineParamsHeader(whichrow, masklist, verbose, csvformat);
5013
5014 if (start < 0) start = 0;
5015 if (count < 0) count = nParam;
5016 int end = start + count;
5017 if (end > nParam) end = nParam;
5018 int paramidoffset = (resetparamid) ? (-start) : 0;
5019
5020 for (int i = start; i < end; ++i) {
5021 if (i > start) {
5022 oss << ",";
5023 }
5024 std::string sFix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
5025 if (csvformat) {
5026 oss << params[i] << sFix;
5027 } else {
5028 oss << " p" << (i+paramidoffset) << sFix << "= " << right << setw(13) << setprecision(6) << params[i];
5029 }
5030 }
5031
5032 if (csvformat) {
5033 oss << ",";
5034 } else {
5035 oss << endl;
5036 }
5037 oss << formatBaselineParamsFooter(rms, nClipped, verbose, csvformat);
5038
5039 return String(oss);
5040 }
5041
5042}
5043
5044std::string Scantable::formatPiecewiseBaselineParams(const std::vector<int>& ranges, const std::vector<float>& params, const std::vector<bool>& fixed, float rms, int nClipped, const std::string& masklist, int whichrow, bool verbose, bool csvformat) const
5045{
5046 int nOutParam = (int)(params.size());
5047 int nPiece = (int)(ranges.size()) - 1;
5048
5049 if (nOutParam < 1) {
5050 return(" Not fitted");
5051 } else if (nPiece < 0) {
5052 return formatBaselineParams(params, fixed, rms, nClipped, masklist, whichrow, verbose, csvformat);
5053 } else if (nPiece < 1) {
5054 return(" Bad count of the piece edge info");
5055 } else if (nOutParam % nPiece != 0) {
5056 return(" Bad count of the output baseline parameters");
5057 } else {
5058
5059 int nParam = nOutParam / nPiece;
5060
5061 ostringstream oss;
5062 oss << formatBaselineParamsHeader(whichrow, masklist, verbose, csvformat);
5063
5064 if (csvformat) {
5065 for (int i = 0; i < nPiece; ++i) {
5066 oss << ranges[i] << "," << (ranges[i+1]-1) << ",";
5067 oss << formatBaselineParams(params, fixed, rms, 0, masklist, whichrow, false, csvformat, i*nParam, nParam, true);
5068 }
5069 } else {
5070 stringstream ss;
5071 ss << ranges[nPiece] << flush;
5072 int wRange = ss.str().size() * 2 + 5;
5073
5074 for (int i = 0; i < nPiece; ++i) {
5075 ss.str("");
5076 ss << " [" << ranges[i] << "," << (ranges[i+1]-1) << "]";
5077 oss << left << setw(wRange) << ss.str();
5078 oss << formatBaselineParams(params, fixed, rms, 0, masklist, whichrow, false, csvformat, i*nParam, nParam, true);
5079 //oss << endl;
5080 }
5081 }
5082
5083 oss << formatBaselineParamsFooter(rms, nClipped, verbose, csvformat);
5084
5085 return String(oss);
5086 }
5087
5088}
5089
5090bool Scantable::hasSameNchanOverIFs()
5091{
5092 int nIF = nif(-1);
5093 int nCh;
5094 int totalPositiveNChan = 0;
5095 int nPositiveNChan = 0;
5096
5097 for (int i = 0; i < nIF; ++i) {
5098 nCh = nchan(i);
5099 if (nCh > 0) {
5100 totalPositiveNChan += nCh;
5101 nPositiveNChan++;
5102 }
5103 }
5104
5105 return (totalPositiveNChan == (nPositiveNChan * nchan(0)));
5106}
5107
5108std::string Scantable::getMaskRangeList(const std::vector<bool>& mask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, bool verbose)
5109{
5110 if (mask.size() <= 0) {
5111 throw(AipsError("The mask elements should be > 0"));
5112 }
5113 int IF = getIF(whichrow);
5114 if (mask.size() != (uInt)nchan(IF)) {
5115 throw(AipsError("Number of channels in scantable != number of mask elements"));
5116 }
5117
5118 if (verbose) {
5119 LogIO logOs(LogOrigin("Scantable", "getMaskRangeList()", WHERE));
5120 logOs << LogIO::WARN << "The current mask window unit is " << coordInfo;
5121 if (!hasSameNchan) {
5122 logOs << endl << "This mask is only valid for IF=" << IF;
5123 }
5124 logOs << LogIO::POST;
5125 }
5126
5127 std::vector<double> abcissa = getAbcissa(whichrow);
5128 std::vector<int> edge = getMaskEdgeIndices(mask);
5129
5130 ostringstream oss;
5131 oss.setf(ios::fixed);
5132 oss << setprecision(1) << "[";
5133 for (uInt i = 0; i < edge.size(); i+=2) {
5134 if (i > 0) oss << ",";
5135 oss << "[" << (float)abcissa[edge[i]] << "," << (float)abcissa[edge[i+1]] << "]";
5136 }
5137 oss << "]" << flush;
5138
5139 return String(oss);
5140}
5141
5142std::vector<int> Scantable::getMaskEdgeIndices(const std::vector<bool>& mask)
5143{
5144 if (mask.size() <= 0) {
5145 throw(AipsError("The mask elements should be > 0"));
5146 }
5147
5148 std::vector<int> out, startIndices, endIndices;
5149 int maskSize = mask.size();
5150
5151 startIndices.clear();
5152 endIndices.clear();
5153
5154 if (mask[0]) {
5155 startIndices.push_back(0);
5156 }
5157 for (int i = 1; i < maskSize; ++i) {
5158 if ((!mask[i-1]) && mask[i]) {
5159 startIndices.push_back(i);
5160 } else if (mask[i-1] && (!mask[i])) {
5161 endIndices.push_back(i-1);
5162 }
5163 }
5164 if (mask[maskSize-1]) {
5165 endIndices.push_back(maskSize-1);
5166 }
5167
5168 if (startIndices.size() != endIndices.size()) {
5169 throw(AipsError("Inconsistent Mask Size: bad data?"));
5170 }
5171 for (uInt i = 0; i < startIndices.size(); ++i) {
5172 if (startIndices[i] > endIndices[i]) {
5173 throw(AipsError("Mask start index > mask end index"));
5174 }
5175 }
5176
5177 out.clear();
5178 for (uInt i = 0; i < startIndices.size(); ++i) {
5179 out.push_back(startIndices[i]);
5180 out.push_back(endIndices[i]);
5181 }
5182
5183 return out;
5184}
5185
5186void Scantable::setTsys(const std::vector<float>& newvals, int whichrow) {
5187 Vector<Float> tsys(newvals);
5188 if (whichrow > -1) {
5189 if (tsysCol_.shape(whichrow) != tsys.shape())
5190 throw(AipsError("Given Tsys values are not of the same shape"));
5191 tsysCol_.put(whichrow, tsys);
5192 } else {
5193 tsysCol_.fillColumn(tsys);
5194 }
5195}
5196
5197vector<float> Scantable::getTsysSpectrum( int whichrow ) const
5198{
5199 Vector<Float> tsys( tsysCol_(whichrow) ) ;
5200 vector<float> stlTsys ;
5201 tsys.tovector( stlTsys ) ;
5202 return stlTsys ;
5203}
5204
5205vector<uint> Scantable::getMoleculeIdColumnData() const
5206{
5207 Vector<uInt> molIds(mmolidCol_.getColumn());
5208 vector<uint> res;
5209 molIds.tovector(res);
5210 return res;
5211}
5212
5213void Scantable::setMoleculeIdColumnData(const std::vector<uint>& molids)
5214{
5215 Vector<uInt> molIds(molids);
5216 Vector<uInt> arr(mmolidCol_.getColumn());
5217 if ( molIds.nelements() != arr.nelements() )
5218 throw AipsError("The input data size must be the number of rows.");
5219 mmolidCol_.putColumn(molIds);
5220}
5221
5222
5223void Scantable::dropXPol()
5224{
5225 if (npol() <= 2) {
5226 return;
5227 }
5228 if (!selector_.empty()) {
5229 throw AipsError("Can only operate with empty selection");
5230 }
5231 std::string taql = "SELECT FROM $1 WHERE POLNO IN [0,1]";
5232 Table tab = tableCommand(taql, table_);
5233 table_ = tab;
5234 table_.rwKeywordSet().define("nPol", Int(2));
5235 originalTable_ = table_;
5236 attach();
5237}
5238
5239}
5240//namespace asap
Note: See TracBrowser for help on using the repository browser.