source: trunk/src/Scantable.cpp@ 2968

Last change on this file since 2968 was 2968, checked in by WataruKawasaki, 11 years ago

New Development: No

JIRA Issue: Yes CAS-6583

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs: test_sdbaseline

Put in Release Notes: No

Module(s): sd

Description: modified *Baseline() to skip baseline fitting/subtraction for row-flagged spectra.


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