source: trunk/src/Scantable.cpp@ 2331

Last change on this file since 2331 was 2321, checked in by Malte Marquarding, 13 years ago

Ticket #249: scantbable schema was changed. This required an update to scantable version 4. I have added a class STUpgrade to handle all schema upgrades. This is handled transparently in the scantable constructor

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 108.6 KB
Line 
1//
2// C++ Implementation: Scantable
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2005
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12#include <map>
13
14#include <atnf/PKSIO/SrcType.h>
15
16#include <casa/aips.h>
17#include <casa/iomanip.h>
18#include <casa/iostream.h>
19#include <casa/OS/File.h>
20#include <casa/OS/Path.h>
21#include <casa/Arrays/Array.h>
22#include <casa/Arrays/ArrayAccessor.h>
23#include <casa/Arrays/ArrayLogical.h>
24#include <casa/Arrays/ArrayMath.h>
25#include <casa/Arrays/MaskArrMath.h>
26#include <casa/Arrays/Slice.h>
27#include <casa/Arrays/Vector.h>
28#include <casa/Arrays/VectorSTLIterator.h>
29#include <casa/BasicMath/Math.h>
30#include <casa/BasicSL/Constants.h>
31#include <casa/Containers/RecordField.h>
32#include <casa/Logging/LogIO.h>
33#include <casa/Quanta/MVAngle.h>
34#include <casa/Quanta/MVTime.h>
35#include <casa/Utilities/GenSort.h>
36
37#include <coordinates/Coordinates/CoordinateUtil.h>
38
39// needed to avoid error in .tcc
40#include <measures/Measures/MCDirection.h>
41//
42#include <measures/Measures/MDirection.h>
43#include <measures/Measures/MEpoch.h>
44#include <measures/Measures/MFrequency.h>
45#include <measures/Measures/MeasRef.h>
46#include <measures/Measures/MeasTable.h>
47#include <measures/TableMeasures/ScalarMeasColumn.h>
48#include <measures/TableMeasures/TableMeasDesc.h>
49#include <measures/TableMeasures/TableMeasRefDesc.h>
50#include <measures/TableMeasures/TableMeasValueDesc.h>
51
52#include <tables/Tables/ArrColDesc.h>
53#include <tables/Tables/ExprNode.h>
54#include <tables/Tables/ScaColDesc.h>
55#include <tables/Tables/SetupNewTab.h>
56#include <tables/Tables/TableCopy.h>
57#include <tables/Tables/TableDesc.h>
58#include <tables/Tables/TableIter.h>
59#include <tables/Tables/TableParse.h>
60#include <tables/Tables/TableRecord.h>
61#include <tables/Tables/TableRow.h>
62#include <tables/Tables/TableVector.h>
63
64#include "MathUtils.h"
65#include "STAttr.h"
66#include "STLineFinder.h"
67#include "STPolCircular.h"
68#include "STPolLinear.h"
69#include "STPolStokes.h"
70#include "STUpgrade.h"
71#include "Scantable.h"
72
73using namespace casa;
74
75namespace asap {
76
77std::map<std::string, STPol::STPolFactory *> Scantable::factories_;
78
79void Scantable::initFactories() {
80 if ( factories_.empty() ) {
81 Scantable::factories_["linear"] = &STPolLinear::myFactory;
82 Scantable::factories_["circular"] = &STPolCircular::myFactory;
83 Scantable::factories_["stokes"] = &STPolStokes::myFactory;
84 }
85}
86
87Scantable::Scantable(Table::TableType ttype) :
88 type_(ttype)
89{
90 initFactories();
91 setupMainTable();
92 freqTable_ = STFrequencies(*this);
93 table_.rwKeywordSet().defineTable("FREQUENCIES", freqTable_.table());
94 weatherTable_ = STWeather(*this);
95 table_.rwKeywordSet().defineTable("WEATHER", weatherTable_.table());
96 focusTable_ = STFocus(*this);
97 table_.rwKeywordSet().defineTable("FOCUS", focusTable_.table());
98 tcalTable_ = STTcal(*this);
99 table_.rwKeywordSet().defineTable("TCAL", tcalTable_.table());
100 moleculeTable_ = STMolecules(*this);
101 table_.rwKeywordSet().defineTable("MOLECULES", moleculeTable_.table());
102 historyTable_ = STHistory(*this);
103 table_.rwKeywordSet().defineTable("HISTORY", historyTable_.table());
104 fitTable_ = STFit(*this);
105 table_.rwKeywordSet().defineTable("FIT", fitTable_.table());
106 table_.tableInfo().setType( "Scantable" ) ;
107 originalTable_ = table_;
108 attach();
109}
110
111Scantable::Scantable(const std::string& name, Table::TableType ttype) :
112 type_(ttype)
113{
114 initFactories();
115
116 Table tab(name, Table::Update);
117 uInt version = tab.keywordSet().asuInt("VERSION");
118 if (version != version_) {
119 STUpgrade upgrader(version_);
120 LogIO os( LogOrigin( "Scantable" ) ) ;
121 os << LogIO::WARN
122 << name << " data format version " << version
123 << " is deprecated" << endl
124 << "Running upgrade."<< endl
125 << LogIO::POST ;
126 std::string outname = upgrader.upgrade(name);
127 os << LogIO::WARN
128 << "Data will be loaded from " << outname << " instead of "
129 << name << LogIO::POST ;
130 tab = Table(outname, Table::Update ) ;
131 }
132 if ( type_ == Table::Memory ) {
133 table_ = tab.copyToMemoryTable(generateName());
134 } else {
135 table_ = tab;
136 }
137 table_.tableInfo().setType( "Scantable" ) ;
138
139 attachSubtables();
140 originalTable_ = table_;
141 attach();
142}
143/*
144Scantable::Scantable(const std::string& name, Table::TableType ttype) :
145 type_(ttype)
146{
147 initFactories();
148 Table tab(name, Table::Update);
149 uInt version = tab.keywordSet().asuInt("VERSION");
150 if (version != version_) {
151 throw(AipsError("Unsupported version of ASAP file."));
152 }
153 if ( type_ == Table::Memory ) {
154 table_ = tab.copyToMemoryTable(generateName());
155 } else {
156 table_ = tab;
157 }
158
159 attachSubtables();
160 originalTable_ = table_;
161 attach();
162}
163*/
164
165Scantable::Scantable( const Scantable& other, bool clear ):
166 Logger()
167{
168 // with or without data
169 String newname = String(generateName());
170 type_ = other.table_.tableType();
171 if ( other.table_.tableType() == Table::Memory ) {
172 if ( clear ) {
173 table_ = TableCopy::makeEmptyMemoryTable(newname,
174 other.table_, True);
175 } else
176 table_ = other.table_.copyToMemoryTable(newname);
177 } else {
178 other.table_.deepCopy(newname, Table::New, False,
179 other.table_.endianFormat(),
180 Bool(clear));
181 table_ = Table(newname, Table::Update);
182 table_.markForDelete();
183 }
184 table_.tableInfo().setType( "Scantable" ) ;
185 /// @todo reindex SCANNO, recompute nbeam, nif, npol
186 if ( clear ) copySubtables(other);
187 attachSubtables();
188 originalTable_ = table_;
189 attach();
190}
191
192void Scantable::copySubtables(const Scantable& other) {
193 Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
194 TableCopy::copyRows(t, other.freqTable_.table());
195 t = table_.rwKeywordSet().asTable("FOCUS");
196 TableCopy::copyRows(t, other.focusTable_.table());
197 t = table_.rwKeywordSet().asTable("WEATHER");
198 TableCopy::copyRows(t, other.weatherTable_.table());
199 t = table_.rwKeywordSet().asTable("TCAL");
200 TableCopy::copyRows(t, other.tcalTable_.table());
201 t = table_.rwKeywordSet().asTable("MOLECULES");
202 TableCopy::copyRows(t, other.moleculeTable_.table());
203 t = table_.rwKeywordSet().asTable("HISTORY");
204 TableCopy::copyRows(t, other.historyTable_.table());
205 t = table_.rwKeywordSet().asTable("FIT");
206 TableCopy::copyRows(t, other.fitTable_.table());
207}
208
209void Scantable::attachSubtables()
210{
211 freqTable_ = STFrequencies(table_);
212 focusTable_ = STFocus(table_);
213 weatherTable_ = STWeather(table_);
214 tcalTable_ = STTcal(table_);
215 moleculeTable_ = STMolecules(table_);
216 historyTable_ = STHistory(table_);
217 fitTable_ = STFit(table_);
218}
219
220Scantable::~Scantable()
221{
222 //cout << "~Scantable() " << this << endl;
223}
224
225void Scantable::setupMainTable()
226{
227 TableDesc td("", "1", TableDesc::Scratch);
228 td.comment() = "An ASAP Scantable";
229 td.rwKeywordSet().define("VERSION", uInt(version_));
230
231 // n Cycles
232 td.addColumn(ScalarColumnDesc<uInt>("SCANNO"));
233 // new index every nBeam x nIF x nPol
234 td.addColumn(ScalarColumnDesc<uInt>("CYCLENO"));
235
236 td.addColumn(ScalarColumnDesc<uInt>("BEAMNO"));
237 td.addColumn(ScalarColumnDesc<uInt>("IFNO"));
238 // linear, circular, stokes
239 td.rwKeywordSet().define("POLTYPE", String("linear"));
240 td.addColumn(ScalarColumnDesc<uInt>("POLNO"));
241
242 td.addColumn(ScalarColumnDesc<uInt>("FREQ_ID"));
243 td.addColumn(ScalarColumnDesc<uInt>("MOLECULE_ID"));
244
245 ScalarColumnDesc<Int> refbeamnoColumn("REFBEAMNO");
246 refbeamnoColumn.setDefault(Int(-1));
247 td.addColumn(refbeamnoColumn);
248
249 ScalarColumnDesc<uInt> flagrowColumn("FLAGROW");
250 flagrowColumn.setDefault(uInt(0));
251 td.addColumn(flagrowColumn);
252
253 td.addColumn(ScalarColumnDesc<Double>("TIME"));
254 TableMeasRefDesc measRef(MEpoch::UTC); // UTC as default
255 TableMeasValueDesc measVal(td, "TIME");
256 TableMeasDesc<MEpoch> mepochCol(measVal, measRef);
257 mepochCol.write(td);
258
259 td.addColumn(ScalarColumnDesc<Double>("INTERVAL"));
260
261 td.addColumn(ScalarColumnDesc<String>("SRCNAME"));
262 // Type of source (on=0, off=1, other=-1)
263 ScalarColumnDesc<Int> stypeColumn("SRCTYPE");
264 stypeColumn.setDefault(Int(-1));
265 td.addColumn(stypeColumn);
266 td.addColumn(ScalarColumnDesc<String>("FIELDNAME"));
267
268 //The actual Data Vectors
269 td.addColumn(ArrayColumnDesc<Float>("SPECTRA"));
270 td.addColumn(ArrayColumnDesc<uChar>("FLAGTRA"));
271 td.addColumn(ArrayColumnDesc<Float>("TSYS"));
272
273 td.addColumn(ArrayColumnDesc<Double>("DIRECTION",
274 IPosition(1,2),
275 ColumnDesc::Direct));
276 TableMeasRefDesc mdirRef(MDirection::J2000); // default
277 TableMeasValueDesc tmvdMDir(td, "DIRECTION");
278 // the TableMeasDesc gives the column a type
279 TableMeasDesc<MDirection> mdirCol(tmvdMDir, mdirRef);
280 // a uder set table type e.g. GALCTIC, B1950 ...
281 td.rwKeywordSet().define("DIRECTIONREF", String("J2000"));
282 // writing create the measure column
283 mdirCol.write(td);
284 td.addColumn(ScalarColumnDesc<Float>("AZIMUTH"));
285 td.addColumn(ScalarColumnDesc<Float>("ELEVATION"));
286 td.addColumn(ScalarColumnDesc<Float>("OPACITY"));
287
288 td.addColumn(ScalarColumnDesc<uInt>("TCAL_ID"));
289 ScalarColumnDesc<Int> fitColumn("FIT_ID");
290 fitColumn.setDefault(Int(-1));
291 td.addColumn(fitColumn);
292
293 td.addColumn(ScalarColumnDesc<uInt>("FOCUS_ID"));
294 td.addColumn(ScalarColumnDesc<uInt>("WEATHER_ID"));
295
296 // columns which just get dragged along, as they aren't used in asap
297 td.addColumn(ScalarColumnDesc<Double>("SRCVELOCITY"));
298 td.addColumn(ArrayColumnDesc<Double>("SRCPROPERMOTION"));
299 td.addColumn(ArrayColumnDesc<Double>("SRCDIRECTION"));
300 td.addColumn(ArrayColumnDesc<Double>("SCANRATE"));
301
302 td.rwKeywordSet().define("OBSMODE", String(""));
303
304 // Now create Table SetUp from the description.
305 SetupNewTable aNewTab(generateName(), td, Table::Scratch);
306 table_ = Table(aNewTab, type_, 0);
307 originalTable_ = table_;
308}
309
310void Scantable::attach()
311{
312 timeCol_.attach(table_, "TIME");
313 srcnCol_.attach(table_, "SRCNAME");
314 srctCol_.attach(table_, "SRCTYPE");
315 specCol_.attach(table_, "SPECTRA");
316 flagsCol_.attach(table_, "FLAGTRA");
317 tsysCol_.attach(table_, "TSYS");
318 cycleCol_.attach(table_,"CYCLENO");
319 scanCol_.attach(table_, "SCANNO");
320 beamCol_.attach(table_, "BEAMNO");
321 ifCol_.attach(table_, "IFNO");
322 polCol_.attach(table_, "POLNO");
323 integrCol_.attach(table_, "INTERVAL");
324 azCol_.attach(table_, "AZIMUTH");
325 elCol_.attach(table_, "ELEVATION");
326 dirCol_.attach(table_, "DIRECTION");
327 fldnCol_.attach(table_, "FIELDNAME");
328 rbeamCol_.attach(table_, "REFBEAMNO");
329
330 mweatheridCol_.attach(table_,"WEATHER_ID");
331 mfitidCol_.attach(table_,"FIT_ID");
332 mfreqidCol_.attach(table_, "FREQ_ID");
333 mtcalidCol_.attach(table_, "TCAL_ID");
334 mfocusidCol_.attach(table_, "FOCUS_ID");
335 mmolidCol_.attach(table_, "MOLECULE_ID");
336
337 //Add auxiliary column for row-based flagging (CAS-1433 Wataru Kawasaki)
338 attachAuxColumnDef(flagrowCol_, "FLAGROW", 0);
339
340}
341
342template<class T, class T2>
343void Scantable::attachAuxColumnDef(ScalarColumn<T>& col,
344 const String& colName,
345 const T2& defValue)
346{
347 try {
348 col.attach(table_, colName);
349 } catch (TableError& err) {
350 String errMesg = err.getMesg();
351 if (errMesg == "Table column " + colName + " is unknown") {
352 table_.addColumn(ScalarColumnDesc<T>(colName));
353 col.attach(table_, colName);
354 col.fillColumn(static_cast<T>(defValue));
355 } else {
356 throw;
357 }
358 } catch (...) {
359 throw;
360 }
361}
362
363template<class T, class T2>
364void Scantable::attachAuxColumnDef(ArrayColumn<T>& col,
365 const String& colName,
366 const Array<T2>& defValue)
367{
368 try {
369 col.attach(table_, colName);
370 } catch (TableError& err) {
371 String errMesg = err.getMesg();
372 if (errMesg == "Table column " + colName + " is unknown") {
373 table_.addColumn(ArrayColumnDesc<T>(colName));
374 col.attach(table_, colName);
375
376 int size = 0;
377 ArrayIterator<T2>& it = defValue.begin();
378 while (it != defValue.end()) {
379 ++size;
380 ++it;
381 }
382 IPosition ip(1, size);
383 Array<T>& arr(ip);
384 for (int i = 0; i < size; ++i)
385 arr[i] = static_cast<T>(defValue[i]);
386
387 col.fillColumn(arr);
388 } else {
389 throw;
390 }
391 } catch (...) {
392 throw;
393 }
394}
395
396void Scantable::setHeader(const STHeader& sdh)
397{
398 table_.rwKeywordSet().define("nIF", sdh.nif);
399 table_.rwKeywordSet().define("nBeam", sdh.nbeam);
400 table_.rwKeywordSet().define("nPol", sdh.npol);
401 table_.rwKeywordSet().define("nChan", sdh.nchan);
402 table_.rwKeywordSet().define("Observer", sdh.observer);
403 table_.rwKeywordSet().define("Project", sdh.project);
404 table_.rwKeywordSet().define("Obstype", sdh.obstype);
405 table_.rwKeywordSet().define("AntennaName", sdh.antennaname);
406 table_.rwKeywordSet().define("AntennaPosition", sdh.antennaposition);
407 table_.rwKeywordSet().define("Equinox", sdh.equinox);
408 table_.rwKeywordSet().define("FreqRefFrame", sdh.freqref);
409 table_.rwKeywordSet().define("FreqRefVal", sdh.reffreq);
410 table_.rwKeywordSet().define("Bandwidth", sdh.bandwidth);
411 table_.rwKeywordSet().define("UTC", sdh.utc);
412 table_.rwKeywordSet().define("FluxUnit", sdh.fluxunit);
413 table_.rwKeywordSet().define("Epoch", sdh.epoch);
414 table_.rwKeywordSet().define("POLTYPE", sdh.poltype);
415}
416
417STHeader Scantable::getHeader() const
418{
419 STHeader sdh;
420 table_.keywordSet().get("nBeam",sdh.nbeam);
421 table_.keywordSet().get("nIF",sdh.nif);
422 table_.keywordSet().get("nPol",sdh.npol);
423 table_.keywordSet().get("nChan",sdh.nchan);
424 table_.keywordSet().get("Observer", sdh.observer);
425 table_.keywordSet().get("Project", sdh.project);
426 table_.keywordSet().get("Obstype", sdh.obstype);
427 table_.keywordSet().get("AntennaName", sdh.antennaname);
428 table_.keywordSet().get("AntennaPosition", sdh.antennaposition);
429 table_.keywordSet().get("Equinox", sdh.equinox);
430 table_.keywordSet().get("FreqRefFrame", sdh.freqref);
431 table_.keywordSet().get("FreqRefVal", sdh.reffreq);
432 table_.keywordSet().get("Bandwidth", sdh.bandwidth);
433 table_.keywordSet().get("UTC", sdh.utc);
434 table_.keywordSet().get("FluxUnit", sdh.fluxunit);
435 table_.keywordSet().get("Epoch", sdh.epoch);
436 table_.keywordSet().get("POLTYPE", sdh.poltype);
437 return sdh;
438}
439
440void Scantable::setSourceType( int stype )
441{
442 if ( stype < 0 || stype > 1 )
443 throw(AipsError("Illegal sourcetype."));
444 TableVector<Int> tabvec(table_, "SRCTYPE");
445 tabvec = Int(stype);
446}
447
448bool Scantable::conformant( const Scantable& other )
449{
450 return this->getHeader().conformant(other.getHeader());
451}
452
453
454
455std::string Scantable::formatSec(Double x) const
456{
457 Double xcop = x;
458 MVTime mvt(xcop/24./3600.); // make days
459
460 if (x < 59.95)
461 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_HM, 7)+"s";
462 else if (x < 3599.95)
463 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_H,7)+" ";
464 else {
465 ostringstream oss;
466 oss << setw(2) << std::right << setprecision(1) << mvt.hour();
467 oss << ":" << mvt.string(MVTime::TIME_CLEAN_NO_H,7) << " ";
468 return String(oss);
469 }
470};
471
472std::string Scantable::formatDirection(const MDirection& md) const
473{
474 Vector<Double> t = md.getAngle(Unit(String("rad"))).getValue();
475 Int prec = 7;
476
477 MVAngle mvLon(t[0]);
478 String sLon = mvLon.string(MVAngle::TIME,prec);
479 uInt tp = md.getRef().getType();
480 if (tp == MDirection::GALACTIC ||
481 tp == MDirection::SUPERGAL ) {
482 sLon = mvLon(0.0).string(MVAngle::ANGLE_CLEAN,prec);
483 }
484 MVAngle mvLat(t[1]);
485 String sLat = mvLat.string(MVAngle::ANGLE+MVAngle::DIG2,prec);
486 return sLon + String(" ") + sLat;
487}
488
489
490std::string Scantable::getFluxUnit() const
491{
492 return table_.keywordSet().asString("FluxUnit");
493}
494
495void Scantable::setFluxUnit(const std::string& unit)
496{
497 String tmp(unit);
498 Unit tU(tmp);
499 if (tU==Unit("K") || tU==Unit("Jy")) {
500 table_.rwKeywordSet().define(String("FluxUnit"), tmp);
501 } else {
502 throw AipsError("Illegal unit - must be compatible with Jy or K");
503 }
504}
505
506void Scantable::setInstrument(const std::string& name)
507{
508 bool throwIt = true;
509 // create an Instrument to see if this is valid
510 STAttr::convertInstrument(name, throwIt);
511 String nameU(name);
512 nameU.upcase();
513 table_.rwKeywordSet().define(String("AntennaName"), nameU);
514}
515
516void Scantable::setFeedType(const std::string& feedtype)
517{
518 if ( Scantable::factories_.find(feedtype) == Scantable::factories_.end() ) {
519 std::string msg = "Illegal feed type "+ feedtype;
520 throw(casa::AipsError(msg));
521 }
522 table_.rwKeywordSet().define(String("POLTYPE"), feedtype);
523}
524
525MPosition Scantable::getAntennaPosition() const
526{
527 Vector<Double> antpos;
528 table_.keywordSet().get("AntennaPosition", antpos);
529 MVPosition mvpos(antpos(0),antpos(1),antpos(2));
530 return MPosition(mvpos);
531}
532
533void Scantable::makePersistent(const std::string& filename)
534{
535 String inname(filename);
536 Path path(inname);
537 /// @todo reindex SCANNO, recompute nbeam, nif, npol
538 inname = path.expandedName();
539 // 2011/03/04 TN
540 // We can comment out this workaround since the essential bug is
541 // fixed in casacore (r20889 in google code).
542 table_.deepCopy(inname, Table::New);
543// // WORKAROUND !!! for Table bug
544// // Remove when fixed in casacore
545// if ( table_.tableType() == Table::Memory && !selector_.empty() ) {
546// Table tab = table_.copyToMemoryTable(generateName());
547// tab.deepCopy(inname, Table::New);
548// tab.markForDelete();
549//
550// } else {
551// table_.deepCopy(inname, Table::New);
552// }
553}
554
555int Scantable::nbeam( int scanno ) const
556{
557 if ( scanno < 0 ) {
558 Int n;
559 table_.keywordSet().get("nBeam",n);
560 return int(n);
561 } else {
562 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
563 Table t = table_(table_.col("SCANNO") == scanno);
564 ROTableRow row(t);
565 const TableRecord& rec = row.get(0);
566 Table subt = t( t.col("IFNO") == Int(rec.asuInt("IFNO"))
567 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
568 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
569 ROTableVector<uInt> v(subt, "BEAMNO");
570 return int(v.nelements());
571 }
572 return 0;
573}
574
575int Scantable::nif( int scanno ) const
576{
577 if ( scanno < 0 ) {
578 Int n;
579 table_.keywordSet().get("nIF",n);
580 return int(n);
581 } else {
582 // take the first POLNO,BEAMNO,CYCLENO as nbeam shouldn't vary with these
583 Table t = table_(table_.col("SCANNO") == scanno);
584 ROTableRow row(t);
585 const TableRecord& rec = row.get(0);
586 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
587 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
588 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
589 if ( subt.nrow() == 0 ) return 0;
590 ROTableVector<uInt> v(subt, "IFNO");
591 return int(v.nelements());
592 }
593 return 0;
594}
595
596int Scantable::npol( int scanno ) const
597{
598 if ( scanno < 0 ) {
599 Int n;
600 table_.keywordSet().get("nPol",n);
601 return n;
602 } else {
603 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
604 Table t = table_(table_.col("SCANNO") == scanno);
605 ROTableRow row(t);
606 const TableRecord& rec = row.get(0);
607 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
608 && t.col("IFNO") == Int(rec.asuInt("IFNO"))
609 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
610 if ( subt.nrow() == 0 ) return 0;
611 ROTableVector<uInt> v(subt, "POLNO");
612 return int(v.nelements());
613 }
614 return 0;
615}
616
617int Scantable::ncycle( int scanno ) const
618{
619 if ( scanno < 0 ) {
620 Block<String> cols(2);
621 cols[0] = "SCANNO";
622 cols[1] = "CYCLENO";
623 TableIterator it(table_, cols);
624 int n = 0;
625 while ( !it.pastEnd() ) {
626 ++n;
627 ++it;
628 }
629 return n;
630 } else {
631 Table t = table_(table_.col("SCANNO") == scanno);
632 ROTableRow row(t);
633 const TableRecord& rec = row.get(0);
634 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
635 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
636 && t.col("IFNO") == Int(rec.asuInt("IFNO")) );
637 if ( subt.nrow() == 0 ) return 0;
638 return int(subt.nrow());
639 }
640 return 0;
641}
642
643
644int Scantable::nrow( int scanno ) const
645{
646 return int(table_.nrow());
647}
648
649int Scantable::nchan( int ifno ) const
650{
651 if ( ifno < 0 ) {
652 Int n;
653 table_.keywordSet().get("nChan",n);
654 return int(n);
655 } else {
656 // take the first SCANNO,POLNO,BEAMNO,CYCLENO as nbeam shouldn't
657 // vary with these
658 Table t = table_(table_.col("IFNO") == ifno, 1);
659 if ( t.nrow() == 0 ) return 0;
660 ROArrayColumn<Float> v(t, "SPECTRA");
661 return v.shape(0)(0);
662 }
663 return 0;
664}
665
666int Scantable::nscan() const {
667 Vector<uInt> scannos(scanCol_.getColumn());
668 uInt nout = genSort( scannos, Sort::Ascending,
669 Sort::QuickSort|Sort::NoDuplicates );
670 return int(nout);
671}
672
673int Scantable::getChannels(int whichrow) const
674{
675 return specCol_.shape(whichrow)(0);
676}
677
678int Scantable::getBeam(int whichrow) const
679{
680 return beamCol_(whichrow);
681}
682
683std::vector<uint> Scantable::getNumbers(const ScalarColumn<uInt>& col) const
684{
685 Vector<uInt> nos(col.getColumn());
686 uInt n = genSort( nos, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates );
687 nos.resize(n, True);
688 std::vector<uint> stlout;
689 nos.tovector(stlout);
690 return stlout;
691}
692
693int Scantable::getIF(int whichrow) const
694{
695 return ifCol_(whichrow);
696}
697
698int Scantable::getPol(int whichrow) const
699{
700 return polCol_(whichrow);
701}
702
703std::string Scantable::formatTime(const MEpoch& me, bool showdate) const
704{
705 return formatTime(me, showdate, 0);
706}
707
708std::string Scantable::formatTime(const MEpoch& me, bool showdate, uInt prec) const
709{
710 MVTime mvt(me.getValue());
711 if (showdate)
712 //mvt.setFormat(MVTime::YMD);
713 mvt.setFormat(MVTime::YMD, prec);
714 else
715 //mvt.setFormat(MVTime::TIME);
716 mvt.setFormat(MVTime::TIME, prec);
717 ostringstream oss;
718 oss << mvt;
719 return String(oss);
720}
721
722void Scantable::calculateAZEL()
723{
724 MPosition mp = getAntennaPosition();
725 MEpoch::ROScalarColumn timeCol(table_, "TIME");
726 ostringstream oss;
727 oss << "Computed azimuth/elevation using " << endl
728 << mp << endl;
729 for (Int i=0; i<nrow(); ++i) {
730 MEpoch me = timeCol(i);
731 MDirection md = getDirection(i);
732 oss << " Time: " << formatTime(me,False) << " Direction: " << formatDirection(md)
733 << endl << " => ";
734 MeasFrame frame(mp, me);
735 Vector<Double> azel =
736 MDirection::Convert(md, MDirection::Ref(MDirection::AZEL,
737 frame)
738 )().getAngle("rad").getValue();
739 azCol_.put(i,Float(azel[0]));
740 elCol_.put(i,Float(azel[1]));
741 oss << "azel: " << azel[0]/C::pi*180.0 << " "
742 << azel[1]/C::pi*180.0 << " (deg)" << endl;
743 }
744 pushLog(String(oss));
745}
746
747void Scantable::clip(const Float uthres, const Float dthres, bool clipoutside, bool unflag)
748{
749 for (uInt i=0; i<table_.nrow(); ++i) {
750 Vector<uChar> flgs = flagsCol_(i);
751 srchChannelsToClip(i, uthres, dthres, clipoutside, unflag, flgs);
752 flagsCol_.put(i, flgs);
753 }
754}
755
756std::vector<bool> Scantable::getClipMask(int whichrow, const Float uthres, const Float dthres, bool clipoutside, bool unflag)
757{
758 Vector<uChar> flags;
759 flagsCol_.get(uInt(whichrow), flags);
760 srchChannelsToClip(uInt(whichrow), uthres, dthres, clipoutside, unflag, flags);
761 Vector<Bool> bflag(flags.shape());
762 convertArray(bflag, flags);
763 //bflag = !bflag;
764
765 std::vector<bool> mask;
766 bflag.tovector(mask);
767 return mask;
768}
769
770void Scantable::srchChannelsToClip(uInt whichrow, const Float uthres, const Float dthres, bool clipoutside, bool unflag,
771 Vector<uChar> flgs)
772{
773 Vector<Float> spcs = specCol_(whichrow);
774 uInt nchannel = nchan();
775 if (spcs.nelements() != nchannel) {
776 throw(AipsError("Data has incorrect number of channels"));
777 }
778 uChar userflag = 1 << 7;
779 if (unflag) {
780 userflag = 0 << 7;
781 }
782 if (clipoutside) {
783 for (uInt j = 0; j < nchannel; ++j) {
784 Float spc = spcs(j);
785 if ((spc >= uthres) || (spc <= dthres)) {
786 flgs(j) = userflag;
787 }
788 }
789 } else {
790 for (uInt j = 0; j < nchannel; ++j) {
791 Float spc = spcs(j);
792 if ((spc < uthres) && (spc > dthres)) {
793 flgs(j) = userflag;
794 }
795 }
796 }
797}
798
799
800void Scantable::flag( int whichrow, const std::vector<bool>& msk, bool unflag ) {
801 std::vector<bool>::const_iterator it;
802 uInt ntrue = 0;
803 if (whichrow >= int(table_.nrow()) ) {
804 throw(AipsError("Invalid row number"));
805 }
806 for (it = msk.begin(); it != msk.end(); ++it) {
807 if ( *it ) {
808 ntrue++;
809 }
810 }
811 //if ( selector_.empty() && (msk.size() == 0 || msk.size() == ntrue) )
812 if ( whichrow == -1 && !unflag && selector_.empty() && (msk.size() == 0 || msk.size() == ntrue) )
813 throw(AipsError("Trying to flag whole scantable."));
814 uChar userflag = 1 << 7;
815 if ( unflag ) {
816 userflag = 0 << 7;
817 }
818 if (whichrow > -1 ) {
819 applyChanFlag(uInt(whichrow), msk, userflag);
820 } else {
821 for ( uInt i=0; i<table_.nrow(); ++i) {
822 applyChanFlag(i, msk, userflag);
823 }
824 }
825}
826
827void Scantable::applyChanFlag( uInt whichrow, const std::vector<bool>& msk, uChar flagval )
828{
829 if (whichrow >= table_.nrow() ) {
830 throw( casa::indexError<int>( whichrow, "asap::Scantable::applyChanFlag: Invalid row number" ) );
831 }
832 Vector<uChar> flgs = flagsCol_(whichrow);
833 if ( msk.size() == 0 ) {
834 flgs = flagval;
835 flagsCol_.put(whichrow, flgs);
836 return;
837 }
838 if ( int(msk.size()) != nchan() ) {
839 throw(AipsError("Mask has incorrect number of channels."));
840 }
841 if ( flgs.nelements() != msk.size() ) {
842 throw(AipsError("Mask has incorrect number of channels."
843 " Probably varying with IF. Please flag per IF"));
844 }
845 std::vector<bool>::const_iterator it;
846 uInt j = 0;
847 for (it = msk.begin(); it != msk.end(); ++it) {
848 if ( *it ) {
849 flgs(j) = flagval;
850 }
851 ++j;
852 }
853 flagsCol_.put(whichrow, flgs);
854}
855
856void Scantable::flagRow(const std::vector<uInt>& rows, bool unflag)
857{
858 if ( selector_.empty() && (rows.size() == table_.nrow()) )
859 throw(AipsError("Trying to flag whole scantable."));
860
861 uInt rowflag = (unflag ? 0 : 1);
862 std::vector<uInt>::const_iterator it;
863 for (it = rows.begin(); it != rows.end(); ++it)
864 flagrowCol_.put(*it, rowflag);
865}
866
867std::vector<bool> Scantable::getMask(int whichrow) const
868{
869 Vector<uChar> flags;
870 flagsCol_.get(uInt(whichrow), flags);
871 Vector<Bool> bflag(flags.shape());
872 convertArray(bflag, flags);
873 bflag = !bflag;
874 std::vector<bool> mask;
875 bflag.tovector(mask);
876 return mask;
877}
878
879std::vector<float> Scantable::getSpectrum( int whichrow,
880 const std::string& poltype ) const
881{
882 String ptype = poltype;
883 if (poltype == "" ) ptype = getPolType();
884 if ( whichrow < 0 || whichrow >= nrow() )
885 throw(AipsError("Illegal row number."));
886 std::vector<float> out;
887 Vector<Float> arr;
888 uInt requestedpol = polCol_(whichrow);
889 String basetype = getPolType();
890 if ( ptype == basetype ) {
891 specCol_.get(whichrow, arr);
892 } else {
893 CountedPtr<STPol> stpol(STPol::getPolClass(Scantable::factories_,
894 basetype));
895 uInt row = uInt(whichrow);
896 stpol->setSpectra(getPolMatrix(row));
897 Float fang,fhand;
898 fang = focusTable_.getTotalAngle(mfocusidCol_(row));
899 fhand = focusTable_.getFeedHand(mfocusidCol_(row));
900 stpol->setPhaseCorrections(fang, fhand);
901 arr = stpol->getSpectrum(requestedpol, ptype);
902 }
903 if ( arr.nelements() == 0 )
904 pushLog("Not enough polarisations present to do the conversion.");
905 arr.tovector(out);
906 return out;
907}
908
909void Scantable::setSpectrum( const std::vector<float>& spec,
910 int whichrow )
911{
912 Vector<Float> spectrum(spec);
913 Vector<Float> arr;
914 specCol_.get(whichrow, arr);
915 if ( spectrum.nelements() != arr.nelements() )
916 throw AipsError("The spectrum has incorrect number of channels.");
917 specCol_.put(whichrow, spectrum);
918}
919
920
921String Scantable::generateName()
922{
923 return (File::newUniqueName("./","temp")).baseName();
924}
925
926const casa::Table& Scantable::table( ) const
927{
928 return table_;
929}
930
931casa::Table& Scantable::table( )
932{
933 return table_;
934}
935
936std::string Scantable::getPolType() const
937{
938 return table_.keywordSet().asString("POLTYPE");
939}
940
941void Scantable::unsetSelection()
942{
943 table_ = originalTable_;
944 attach();
945 selector_.reset();
946}
947
948void Scantable::setSelection( const STSelector& selection )
949{
950 Table tab = const_cast<STSelector&>(selection).apply(originalTable_);
951 if ( tab.nrow() == 0 ) {
952 throw(AipsError("Selection contains no data. Not applying it."));
953 }
954 table_ = tab;
955 attach();
956// tab.rwKeywordSet().define("nBeam",(Int)(getBeamNos().size())) ;
957// vector<uint> selectedIFs = getIFNos() ;
958// Int newnIF = selectedIFs.size() ;
959// tab.rwKeywordSet().define("nIF",newnIF) ;
960// if ( newnIF != 0 ) {
961// Int newnChan = 0 ;
962// for ( Int i = 0 ; i < newnIF ; i++ ) {
963// Int nChan = nchan( selectedIFs[i] ) ;
964// if ( newnChan > nChan )
965// newnChan = nChan ;
966// }
967// tab.rwKeywordSet().define("nChan",newnChan) ;
968// }
969// tab.rwKeywordSet().define("nPol",(Int)(getPolNos().size())) ;
970 selector_ = selection;
971}
972
973
974std::string Scantable::headerSummary()
975{
976 // Format header info
977// STHeader sdh;
978// sdh = getHeader();
979// sdh.print();
980 ostringstream oss;
981 oss.flags(std::ios_base::left);
982 String tmp;
983 // Project
984 table_.keywordSet().get("Project", tmp);
985 oss << setw(15) << "Project:" << tmp << endl;
986 // Observation date
987 oss << setw(15) << "Obs Date:" << getTime(-1,true) << endl;
988 // Observer
989 oss << setw(15) << "Observer:"
990 << table_.keywordSet().asString("Observer") << endl;
991 // Antenna Name
992 table_.keywordSet().get("AntennaName", tmp);
993 oss << setw(15) << "Antenna Name:" << tmp << endl;
994 // Obs type
995 table_.keywordSet().get("Obstype", tmp);
996 // Records (nrow)
997 oss << setw(15) << "Data Records:" << table_.nrow() << " rows" << endl;
998 oss << setw(15) << "Obs. Type:" << tmp << endl;
999 // Beams, IFs, Polarizations, and Channels
1000 oss << setw(15) << "Beams:" << setw(4) << nbeam() << endl
1001 << setw(15) << "IFs:" << setw(4) << nif() << endl
1002 << setw(15) << "Polarisations:" << setw(4) << npol()
1003 << "(" << getPolType() << ")" << endl
1004 << setw(15) << "Channels:" << nchan() << endl;
1005 // Flux unit
1006 table_.keywordSet().get("FluxUnit", tmp);
1007 oss << setw(15) << "Flux Unit:" << tmp << endl;
1008 // Abscissa Unit
1009 oss << setw(15) << "Abscissa:" << getAbcissaLabel(0) << endl;
1010 // Selection
1011 oss << selector_.print() << endl;
1012
1013 return String(oss);
1014}
1015
1016void Scantable::summary( const std::string& filename )
1017{
1018 ostringstream oss;
1019 ofstream ofs;
1020 LogIO ols(LogOrigin("Scantable", "summary", WHERE));
1021
1022 if (filename != "")
1023 ofs.open( filename.c_str(), ios::out );
1024
1025 oss << endl;
1026 oss << asap::SEPERATOR << endl;
1027 oss << " Scan Table Summary" << endl;
1028 oss << asap::SEPERATOR << endl;
1029
1030 // Format header info
1031 oss << headerSummary();
1032 oss << endl;
1033
1034 if (table_.nrow() <= 0){
1035 oss << asap::SEPERATOR << endl;
1036 oss << "The MAIN table is empty: there are no data!!!" << endl;
1037 oss << asap::SEPERATOR << endl;
1038
1039 ols << String(oss) << LogIO::POST;
1040 if (ofs) {
1041 ofs << String(oss) << flush;
1042 ofs.close();
1043 }
1044 return;
1045 }
1046
1047
1048
1049 // main table
1050 String dirtype = "Position ("
1051 + getDirectionRefString()
1052 + ")";
1053 oss.flags(std::ios_base::left);
1054 oss << setw(5) << "Scan"
1055 << setw(15) << "Source"
1056 << setw(35) << "Time range"
1057 << setw(2) << "" << setw(7) << "Int[s]"
1058 << setw(7) << "Record"
1059 << setw(8) << "SrcType"
1060 << setw(8) << "FreqIDs"
1061 << setw(7) << "MolIDs" << endl;
1062 oss << setw(7)<< "" << setw(6) << "Beam"
1063 << setw(23) << dirtype << endl;
1064
1065 oss << asap::SEPERATOR << endl;
1066
1067 // Flush summary and clear up the string
1068 ols << String(oss) << LogIO::POST;
1069 if (ofs) ofs << String(oss) << flush;
1070 oss.str("");
1071 oss.clear();
1072
1073
1074 // Get Freq_ID map
1075 ROScalarColumn<uInt> ftabIds(frequencies().table(), "ID");
1076 Int nfid = ftabIds.nrow();
1077 if (nfid <= 0){
1078 oss << "FREQUENCIES subtable is empty: there are no data!!!" << endl;
1079 oss << asap::SEPERATOR << endl;
1080
1081 ols << String(oss) << LogIO::POST;
1082 if (ofs) {
1083 ofs << String(oss) << flush;
1084 ofs.close();
1085 }
1086 return;
1087 }
1088 // Storages of overall IFNO, POLNO, and nchan per FREQ_ID
1089 // the orders are identical to ID in FREQ subtable
1090 Block< Vector<uInt> > ifNos(nfid), polNos(nfid);
1091 Vector<Int> fIdchans(nfid,-1);
1092 map<uInt, Int> fidMap; // (FREQ_ID, row # in FREQ subtable) pair
1093 for (Int i=0; i < nfid; i++){
1094 // fidMap[freqId] returns row number in FREQ subtable
1095 fidMap.insert(pair<uInt, Int>(ftabIds(i),i));
1096 ifNos[i] = Vector<uInt>();
1097 polNos[i] = Vector<uInt>();
1098 }
1099
1100 TableIterator iter(table_, "SCANNO");
1101
1102 // Vars for keeping track of time, freqids, molIds in a SCANNO
1103 Vector<uInt> freqids;
1104 Vector<uInt> molids;
1105 Vector<uInt> beamids(1,0);
1106 Vector<MDirection> beamDirs;
1107 Vector<Int> stypeids(1,0);
1108 Vector<String> stypestrs;
1109 Int nfreq(1);
1110 Int nmol(1);
1111 uInt nbeam(1);
1112 uInt nstype(1);
1113
1114 Double btime(0.0), etime(0.0);
1115 Double meanIntTim(0.0);
1116
1117 uInt currFreqId(0), ftabRow(0);
1118 Int iflen(0), pollen(0);
1119
1120 while (!iter.pastEnd()) {
1121 Table subt = iter.table();
1122 uInt snrow = subt.nrow();
1123 ROTableRow row(subt);
1124 const TableRecord& rec = row.get(0);
1125
1126 // relevant columns
1127 ROScalarColumn<Double> mjdCol(subt,"TIME");
1128 ROScalarColumn<Double> intervalCol(subt,"INTERVAL");
1129 MDirection::ROScalarColumn dirCol(subt,"DIRECTION");
1130
1131 ScalarColumn<uInt> freqIdCol(subt,"FREQ_ID");
1132 ScalarColumn<uInt> molIdCol(subt,"MOLECULE_ID");
1133 ROScalarColumn<uInt> beamCol(subt,"BEAMNO");
1134 ROScalarColumn<Int> stypeCol(subt,"SRCTYPE");
1135
1136 ROScalarColumn<uInt> ifNoCol(subt,"IFNO");
1137 ROScalarColumn<uInt> polNoCol(subt,"POLNO");
1138
1139
1140 // Times
1141 meanIntTim = sum(intervalCol.getColumn()) / (double) snrow;
1142 minMax(btime, etime, mjdCol.getColumn());
1143 etime += meanIntTim/C::day;
1144
1145 // MOLECULE_ID and FREQ_ID
1146 molids = getNumbers(molIdCol);
1147 molids.shape(nmol);
1148
1149 freqids = getNumbers(freqIdCol);
1150 freqids.shape(nfreq);
1151
1152 // Add first beamid, and srcNames
1153 beamids.resize(1,False);
1154 beamDirs.resize(1,False);
1155 beamids(0)=beamCol(0);
1156 beamDirs(0)=dirCol(0);
1157 nbeam = 1;
1158
1159 stypeids.resize(1,False);
1160 stypeids(0)=stypeCol(0);
1161 nstype = 1;
1162
1163 // Global listings of nchan/IFNO/POLNO per FREQ_ID
1164 currFreqId=freqIdCol(0);
1165 ftabRow = fidMap[currFreqId];
1166 // Assumes an identical number of channels per FREQ_ID
1167 if (fIdchans(ftabRow) < 0 ) {
1168 RORecordFieldPtr< Array<Float> > spec(rec, "SPECTRA");
1169 fIdchans(ftabRow)=(*spec).shape()(0);
1170 }
1171 // Should keep ifNos and polNos form the previous SCANNO
1172 if ( !anyEQ(ifNos[ftabRow],ifNoCol(0)) ) {
1173 ifNos[ftabRow].shape(iflen);
1174 iflen++;
1175 ifNos[ftabRow].resize(iflen,True);
1176 ifNos[ftabRow](iflen-1) = ifNoCol(0);
1177 }
1178 if ( !anyEQ(polNos[ftabRow],polNoCol(0)) ) {
1179 polNos[ftabRow].shape(pollen);
1180 pollen++;
1181 polNos[ftabRow].resize(pollen,True);
1182 polNos[ftabRow](pollen-1) = polNoCol(0);
1183 }
1184
1185 for (uInt i=1; i < snrow; i++){
1186 // Need to list BEAMNO and DIRECTION in the same order
1187 if ( !anyEQ(beamids,beamCol(i)) ) {
1188 nbeam++;
1189 beamids.resize(nbeam,True);
1190 beamids(nbeam-1)=beamCol(i);
1191 beamDirs.resize(nbeam,True);
1192 beamDirs(nbeam-1)=dirCol(i);
1193 }
1194
1195 // SRCTYPE is Int (getNumber takes only uInt)
1196 if ( !anyEQ(stypeids,stypeCol(i)) ) {
1197 nstype++;
1198 stypeids.resize(nstype,True);
1199 stypeids(nstype-1)=stypeCol(i);
1200 }
1201
1202 // Global listings of nchan/IFNO/POLNO per FREQ_ID
1203 currFreqId=freqIdCol(i);
1204 ftabRow = fidMap[currFreqId];
1205 if (fIdchans(ftabRow) < 0 ) {
1206 const TableRecord& rec = row.get(i);
1207 RORecordFieldPtr< Array<Float> > spec(rec, "SPECTRA");
1208 fIdchans(ftabRow) = (*spec).shape()(0);
1209 }
1210 if ( !anyEQ(ifNos[ftabRow],ifNoCol(i)) ) {
1211 ifNos[ftabRow].shape(iflen);
1212 iflen++;
1213 ifNos[ftabRow].resize(iflen,True);
1214 ifNos[ftabRow](iflen-1) = ifNoCol(i);
1215 }
1216 if ( !anyEQ(polNos[ftabRow],polNoCol(i)) ) {
1217 polNos[ftabRow].shape(pollen);
1218 pollen++;
1219 polNos[ftabRow].resize(pollen,True);
1220 polNos[ftabRow](pollen-1) = polNoCol(i);
1221 }
1222 } // end of row iteration
1223
1224 stypestrs.resize(nstype,False);
1225 for (uInt j=0; j < nstype; j++)
1226 stypestrs(j) = SrcType::getName(stypeids(j));
1227
1228 // Format Scan summary
1229 oss << setw(4) << std::right << rec.asuInt("SCANNO")
1230 << std::left << setw(1) << ""
1231 << setw(15) << rec.asString("SRCNAME")
1232 << setw(21) << MVTime(btime).string(MVTime::YMD,7)
1233 << setw(3) << " - " << MVTime(etime).string(MVTime::TIME,7)
1234 << setw(3) << "" << setw(6) << meanIntTim << setw(1) << ""
1235 << std::right << setw(5) << snrow << setw(2) << ""
1236 << std::left << stypestrs << setw(1) << ""
1237 << freqids << setw(1) << ""
1238 << molids << endl;
1239 // Format Beam summary
1240 for (uInt j=0; j < nbeam; j++) {
1241 oss << setw(7) << "" << setw(6) << beamids(j) << setw(1) << ""
1242 << formatDirection(beamDirs(j)) << endl;
1243 }
1244 // Flush summary every scan and clear up the string
1245 ols << String(oss) << LogIO::POST;
1246 if (ofs) ofs << String(oss) << flush;
1247 oss.str("");
1248 oss.clear();
1249
1250 ++iter;
1251 } // end of scan iteration
1252 oss << asap::SEPERATOR << endl;
1253
1254 // List FRECUENCIES Table (using STFrequencies.print may be slow)
1255 oss << "FREQUENCIES: " << nfreq << endl;
1256 oss << std::right << setw(5) << "ID" << setw(2) << ""
1257 << std::left << setw(5) << "IFNO" << setw(2) << ""
1258 << setw(8) << "Frame"
1259 << setw(16) << "RefVal"
1260 << setw(7) << "RefPix"
1261 << setw(15) << "Increment"
1262 << setw(9) << "Channels"
1263 << setw(6) << "POLNOs" << endl;
1264 Int tmplen;
1265 for (Int i=0; i < nfid; i++){
1266 // List row=i of FREQUENCIES subtable
1267 ifNos[i].shape(tmplen);
1268 if (tmplen == 1) {
1269 oss << std::right << setw(5) << ftabIds(i) << setw(2) << ""
1270 << setw(3) << ifNos[i](0) << setw(1) << ""
1271 << std::left << setw(46) << frequencies().print(ftabIds(i))
1272 << setw(2) << ""
1273 << std::right << setw(8) << fIdchans[i] << setw(2) << ""
1274 << std::left << polNos[i] << endl;
1275 } else if (tmplen > 0 ) {
1276 // You shouldn't come here
1277 oss << std::left
1278 << "Multiple IFNOs in FREQ_ID = " << ftabIds(i)
1279 << " !!!" << endl;
1280 }
1281 }
1282 oss << asap::SEPERATOR << endl;
1283
1284 // List MOLECULES Table (currently lists all rows)
1285 oss << "MOLECULES: " << endl;
1286 if (molecules().nrow() <= 0) {
1287 oss << " MOLECULES subtable is empty: there are no data" << endl;
1288 } else {
1289 ROTableRow row(molecules().table());
1290 oss << std::right << setw(5) << "ID"
1291 << std::left << setw(3) << ""
1292 << setw(18) << "RestFreq"
1293 << setw(15) << "Name" << endl;
1294 for (Int i=0; i < molecules().nrow(); i++){
1295 const TableRecord& rec=row.get(i);
1296 oss << std::right << setw(5) << rec.asuInt("ID")
1297 << std::left << setw(3) << ""
1298 << rec.asArrayDouble("RESTFREQUENCY") << setw(1) << ""
1299 << rec.asArrayString("NAME") << endl;
1300 }
1301 }
1302 oss << asap::SEPERATOR << endl;
1303 ols << String(oss) << LogIO::POST;
1304 if (ofs) {
1305 ofs << String(oss) << flush;
1306 ofs.close();
1307 }
1308 // return String(oss);
1309}
1310
1311
1312std::string Scantable::oldheaderSummary()
1313{
1314 // Format header info
1315// STHeader sdh;
1316// sdh = getHeader();
1317// sdh.print();
1318 ostringstream oss;
1319 oss.flags(std::ios_base::left);
1320 oss << setw(15) << "Beams:" << setw(4) << nbeam() << endl
1321 << setw(15) << "IFs:" << setw(4) << nif() << endl
1322 << setw(15) << "Polarisations:" << setw(4) << npol()
1323 << "(" << getPolType() << ")" << endl
1324 << setw(15) << "Channels:" << nchan() << endl;
1325 String tmp;
1326 oss << setw(15) << "Observer:"
1327 << table_.keywordSet().asString("Observer") << endl;
1328 oss << setw(15) << "Obs Date:" << getTime(-1,true) << endl;
1329 table_.keywordSet().get("Project", tmp);
1330 oss << setw(15) << "Project:" << tmp << endl;
1331 table_.keywordSet().get("Obstype", tmp);
1332 oss << setw(15) << "Obs. Type:" << tmp << endl;
1333 table_.keywordSet().get("AntennaName", tmp);
1334 oss << setw(15) << "Antenna Name:" << tmp << endl;
1335 table_.keywordSet().get("FluxUnit", tmp);
1336 oss << setw(15) << "Flux Unit:" << tmp << endl;
1337 int nid = moleculeTable_.nrow();
1338 Bool firstline = True;
1339 oss << setw(15) << "Rest Freqs:";
1340 for (int i=0; i<nid; i++) {
1341 Table t = table_(table_.col("MOLECULE_ID") == i, 1);
1342 if (t.nrow() > 0) {
1343 Vector<Double> vec(moleculeTable_.getRestFrequency(i));
1344 if (vec.nelements() > 0) {
1345 if (firstline) {
1346 oss << setprecision(10) << vec << " [Hz]" << endl;
1347 firstline=False;
1348 }
1349 else{
1350 oss << setw(15)<<" " << setprecision(10) << vec << " [Hz]" << endl;
1351 }
1352 } else {
1353 oss << "none" << endl;
1354 }
1355 }
1356 }
1357
1358 oss << setw(15) << "Abcissa:" << getAbcissaLabel(0) << endl;
1359 oss << selector_.print() << endl;
1360 return String(oss);
1361}
1362
1363 //std::string Scantable::summary( const std::string& filename )
1364void Scantable::oldsummary( const std::string& filename )
1365{
1366 ostringstream oss;
1367 ofstream ofs;
1368 LogIO ols(LogOrigin("Scantable", "summary", WHERE));
1369
1370 if (filename != "")
1371 ofs.open( filename.c_str(), ios::out );
1372
1373 oss << endl;
1374 oss << asap::SEPERATOR << endl;
1375 oss << " Scan Table Summary" << endl;
1376 oss << asap::SEPERATOR << endl;
1377
1378 // Format header info
1379 oss << oldheaderSummary();
1380 oss << endl;
1381
1382 // main table
1383 String dirtype = "Position ("
1384 + getDirectionRefString()
1385 + ")";
1386 oss.flags(std::ios_base::left);
1387 oss << setw(5) << "Scan" << setw(15) << "Source"
1388 << setw(10) << "Time" << setw(18) << "Integration"
1389 << setw(15) << "Source Type" << endl;
1390 oss << setw(5) << "" << setw(5) << "Beam" << setw(3) << "" << dirtype << endl;
1391 oss << setw(10) << "" << setw(3) << "IF" << setw(3) << ""
1392 << setw(8) << "Frame" << setw(16)
1393 << "RefVal" << setw(10) << "RefPix" << setw(12) << "Increment"
1394 << setw(7) << "Channels"
1395 << endl;
1396 oss << asap::SEPERATOR << endl;
1397
1398 // Flush summary and clear up the string
1399 ols << String(oss) << LogIO::POST;
1400 if (ofs) ofs << String(oss) << flush;
1401 oss.str("");
1402 oss.clear();
1403
1404 TableIterator iter(table_, "SCANNO");
1405 while (!iter.pastEnd()) {
1406 Table subt = iter.table();
1407 ROTableRow row(subt);
1408 MEpoch::ROScalarColumn timeCol(subt,"TIME");
1409 const TableRecord& rec = row.get(0);
1410 oss << setw(4) << std::right << rec.asuInt("SCANNO")
1411 << std::left << setw(1) << ""
1412 << setw(15) << rec.asString("SRCNAME")
1413 << setw(10) << formatTime(timeCol(0), false);
1414 // count the cycles in the scan
1415 TableIterator cyciter(subt, "CYCLENO");
1416 int nint = 0;
1417 while (!cyciter.pastEnd()) {
1418 ++nint;
1419 ++cyciter;
1420 }
1421 oss << setw(3) << std::right << nint << setw(3) << " x " << std::left
1422 << setw(11) << formatSec(rec.asFloat("INTERVAL")) << setw(1) << ""
1423 << setw(15) << SrcType::getName(rec.asInt("SRCTYPE")) << endl;
1424
1425 TableIterator biter(subt, "BEAMNO");
1426 while (!biter.pastEnd()) {
1427 Table bsubt = biter.table();
1428 ROTableRow brow(bsubt);
1429 const TableRecord& brec = brow.get(0);
1430 uInt row0 = bsubt.rowNumbers(table_)[0];
1431 oss << setw(5) << "" << setw(4) << std::right << brec.asuInt("BEAMNO")<< std::left;
1432 oss << setw(4) << "" << formatDirection(getDirection(row0)) << endl;
1433 TableIterator iiter(bsubt, "IFNO");
1434 while (!iiter.pastEnd()) {
1435 Table isubt = iiter.table();
1436 ROTableRow irow(isubt);
1437 const TableRecord& irec = irow.get(0);
1438 oss << setw(9) << "";
1439 oss << setw(3) << std::right << irec.asuInt("IFNO") << std::left
1440 << setw(1) << "" << frequencies().print(irec.asuInt("FREQ_ID"))
1441 << setw(3) << "" << nchan(irec.asuInt("IFNO"))
1442 << endl;
1443
1444 ++iiter;
1445 }
1446 ++biter;
1447 }
1448 // Flush summary every scan and clear up the string
1449 ols << String(oss) << LogIO::POST;
1450 if (ofs) ofs << String(oss) << flush;
1451 oss.str("");
1452 oss.clear();
1453
1454 ++iter;
1455 }
1456 oss << asap::SEPERATOR << endl;
1457 ols << String(oss) << LogIO::POST;
1458 if (ofs) {
1459 ofs << String(oss) << flush;
1460 ofs.close();
1461 }
1462 // return String(oss);
1463}
1464
1465// std::string Scantable::getTime(int whichrow, bool showdate) const
1466// {
1467// MEpoch::ROScalarColumn timeCol(table_, "TIME");
1468// MEpoch me;
1469// if (whichrow > -1) {
1470// me = timeCol(uInt(whichrow));
1471// } else {
1472// Double tm;
1473// table_.keywordSet().get("UTC",tm);
1474// me = MEpoch(MVEpoch(tm));
1475// }
1476// return formatTime(me, showdate);
1477// }
1478
1479std::string Scantable::getTime(int whichrow, bool showdate, uInt prec) const
1480{
1481 MEpoch me;
1482 me = getEpoch(whichrow);
1483 return formatTime(me, showdate, prec);
1484}
1485
1486MEpoch Scantable::getEpoch(int whichrow) const
1487{
1488 if (whichrow > -1) {
1489 return timeCol_(uInt(whichrow));
1490 } else {
1491 Double tm;
1492 table_.keywordSet().get("UTC",tm);
1493 return MEpoch(MVEpoch(tm));
1494 }
1495}
1496
1497std::string Scantable::getDirectionString(int whichrow) const
1498{
1499 return formatDirection(getDirection(uInt(whichrow)));
1500}
1501
1502
1503SpectralCoordinate Scantable::getSpectralCoordinate(int whichrow) const {
1504 const MPosition& mp = getAntennaPosition();
1505 const MDirection& md = getDirection(whichrow);
1506 const MEpoch& me = timeCol_(whichrow);
1507 //Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
1508 Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
1509 return freqTable_.getSpectralCoordinate(md, mp, me, rf,
1510 mfreqidCol_(whichrow));
1511}
1512
1513std::vector< double > Scantable::getAbcissa( int whichrow ) const
1514{
1515 if ( whichrow > int(table_.nrow()) ) throw(AipsError("Illegal row number"));
1516 std::vector<double> stlout;
1517 int nchan = specCol_(whichrow).nelements();
1518 String us = freqTable_.getUnitString();
1519 if ( us == "" || us == "pixel" || us == "channel" ) {
1520 for (int i=0; i<nchan; ++i) {
1521 stlout.push_back(double(i));
1522 }
1523 return stlout;
1524 }
1525 SpectralCoordinate spc = getSpectralCoordinate(whichrow);
1526 Vector<Double> pixel(nchan);
1527 Vector<Double> world;
1528 indgen(pixel);
1529 if ( Unit(us) == Unit("Hz") ) {
1530 for ( int i=0; i < nchan; ++i) {
1531 Double world;
1532 spc.toWorld(world, pixel[i]);
1533 stlout.push_back(double(world));
1534 }
1535 } else if ( Unit(us) == Unit("km/s") ) {
1536 Vector<Double> world;
1537 spc.pixelToVelocity(world, pixel);
1538 world.tovector(stlout);
1539 }
1540 return stlout;
1541}
1542void Scantable::setDirectionRefString( const std::string & refstr )
1543{
1544 MDirection::Types mdt;
1545 if (refstr != "" && !MDirection::getType(mdt, refstr)) {
1546 throw(AipsError("Illegal Direction frame."));
1547 }
1548 if ( refstr == "" ) {
1549 String defaultstr = MDirection::showType(dirCol_.getMeasRef().getType());
1550 table_.rwKeywordSet().define("DIRECTIONREF", defaultstr);
1551 } else {
1552 table_.rwKeywordSet().define("DIRECTIONREF", String(refstr));
1553 }
1554}
1555
1556std::string Scantable::getDirectionRefString( ) const
1557{
1558 return table_.keywordSet().asString("DIRECTIONREF");
1559}
1560
1561MDirection Scantable::getDirection(int whichrow ) const
1562{
1563 String usertype = table_.keywordSet().asString("DIRECTIONREF");
1564 String type = MDirection::showType(dirCol_.getMeasRef().getType());
1565 if ( usertype != type ) {
1566 MDirection::Types mdt;
1567 if (!MDirection::getType(mdt, usertype)) {
1568 throw(AipsError("Illegal Direction frame."));
1569 }
1570 return dirCol_.convert(uInt(whichrow), mdt);
1571 } else {
1572 return dirCol_(uInt(whichrow));
1573 }
1574}
1575
1576std::string Scantable::getAbcissaLabel( int whichrow ) const
1577{
1578 if ( whichrow > int(table_.nrow()) ) throw(AipsError("Illegal ro number"));
1579 const MPosition& mp = getAntennaPosition();
1580 const MDirection& md = getDirection(whichrow);
1581 const MEpoch& me = timeCol_(whichrow);
1582 //const Double& rf = mmolidCol_(whichrow);
1583 const Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
1584 SpectralCoordinate spc =
1585 freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
1586
1587 String s = "Channel";
1588 Unit u = Unit(freqTable_.getUnitString());
1589 if (u == Unit("km/s")) {
1590 s = CoordinateUtil::axisLabel(spc, 0, True,True, True);
1591 } else if (u == Unit("Hz")) {
1592 Vector<String> wau(1);wau = u.getName();
1593 spc.setWorldAxisUnits(wau);
1594 s = CoordinateUtil::axisLabel(spc, 0, True, True, False);
1595 }
1596 return s;
1597
1598}
1599
1600/**
1601void asap::Scantable::setRestFrequencies( double rf, const std::string& name,
1602 const std::string& unit )
1603**/
1604void Scantable::setRestFrequencies( vector<double> rf, const vector<std::string>& name,
1605 const std::string& unit )
1606
1607{
1608 ///@todo lookup in line table to fill in name and formattedname
1609 Unit u(unit);
1610 //Quantum<Double> urf(rf, u);
1611 Quantum<Vector<Double> >urf(rf, u);
1612 Vector<String> formattedname(0);
1613 //cerr<<"Scantable::setRestFrequnecies="<<urf<<endl;
1614
1615 //uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), name, "");
1616 uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), mathutil::toVectorString(name), formattedname);
1617 TableVector<uInt> tabvec(table_, "MOLECULE_ID");
1618 tabvec = id;
1619}
1620
1621/**
1622void asap::Scantable::setRestFrequencies( const std::string& name )
1623{
1624 throw(AipsError("setRestFrequencies( const std::string& name ) NYI"));
1625 ///@todo implement
1626}
1627**/
1628
1629void Scantable::setRestFrequencies( const vector<std::string>& name )
1630{
1631 (void) name; // suppress unused warning
1632 throw(AipsError("setRestFrequencies( const vector<std::string>& name ) NYI"));
1633 ///@todo implement
1634}
1635
1636std::vector< unsigned int > Scantable::rownumbers( ) const
1637{
1638 std::vector<unsigned int> stlout;
1639 Vector<uInt> vec = table_.rowNumbers();
1640 vec.tovector(stlout);
1641 return stlout;
1642}
1643
1644
1645Matrix<Float> Scantable::getPolMatrix( uInt whichrow ) const
1646{
1647 ROTableRow row(table_);
1648 const TableRecord& rec = row.get(whichrow);
1649 Table t =
1650 originalTable_( originalTable_.col("SCANNO") == Int(rec.asuInt("SCANNO"))
1651 && originalTable_.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
1652 && originalTable_.col("IFNO") == Int(rec.asuInt("IFNO"))
1653 && originalTable_.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
1654 ROArrayColumn<Float> speccol(t, "SPECTRA");
1655 return speccol.getColumn();
1656}
1657
1658std::vector< std::string > Scantable::columnNames( ) const
1659{
1660 Vector<String> vec = table_.tableDesc().columnNames();
1661 return mathutil::tovectorstring(vec);
1662}
1663
1664MEpoch::Types Scantable::getTimeReference( ) const
1665{
1666 return MEpoch::castType(timeCol_.getMeasRef().getType());
1667}
1668
1669void Scantable::addFit( const STFitEntry& fit, int row )
1670{
1671 //cout << mfitidCol_(uInt(row)) << endl;
1672 LogIO os( LogOrigin( "Scantable", "addFit()", WHERE ) ) ;
1673 os << mfitidCol_(uInt(row)) << LogIO::POST ;
1674 uInt id = fitTable_.addEntry(fit, mfitidCol_(uInt(row)));
1675 mfitidCol_.put(uInt(row), id);
1676}
1677
1678void Scantable::shift(int npix)
1679{
1680 Vector<uInt> fids(mfreqidCol_.getColumn());
1681 genSort( fids, Sort::Ascending,
1682 Sort::QuickSort|Sort::NoDuplicates );
1683 for (uInt i=0; i<fids.nelements(); ++i) {
1684 frequencies().shiftRefPix(npix, fids[i]);
1685 }
1686}
1687
1688String Scantable::getAntennaName() const
1689{
1690 String out;
1691 table_.keywordSet().get("AntennaName", out);
1692 String::size_type pos1 = out.find("@") ;
1693 String::size_type pos2 = out.find("//") ;
1694 if ( pos2 != String::npos )
1695 out = out.substr(pos2+2,pos1-pos2-2) ;
1696 else if ( pos1 != String::npos )
1697 out = out.substr(0,pos1) ;
1698 return out;
1699}
1700
1701int Scantable::checkScanInfo(const std::vector<int>& scanlist) const
1702{
1703 String tbpath;
1704 int ret = 0;
1705 if ( table_.keywordSet().isDefined("GBT_GO") ) {
1706 table_.keywordSet().get("GBT_GO", tbpath);
1707 Table t(tbpath,Table::Old);
1708 // check each scan if other scan of the pair exist
1709 int nscan = scanlist.size();
1710 for (int i = 0; i < nscan; i++) {
1711 Table subt = t( t.col("SCAN") == scanlist[i]+1 );
1712 if (subt.nrow()==0) {
1713 //cerr <<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<endl;
1714 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1715 os <<LogIO::WARN<<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<LogIO::POST;
1716 ret = 1;
1717 break;
1718 }
1719 ROTableRow row(subt);
1720 const TableRecord& rec = row.get(0);
1721 int scan1seqn = rec.asuInt("PROCSEQN");
1722 int laston1 = rec.asuInt("LASTON");
1723 if ( rec.asuInt("PROCSIZE")==2 ) {
1724 if ( i < nscan-1 ) {
1725 Table subt2 = t( t.col("SCAN") == scanlist[i+1]+1 );
1726 if ( subt2.nrow() == 0) {
1727 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1728
1729 //cerr<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<endl;
1730 os<<LogIO::WARN<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<LogIO::POST;
1731 ret = 1;
1732 break;
1733 }
1734 ROTableRow row2(subt2);
1735 const TableRecord& rec2 = row2.get(0);
1736 int scan2seqn = rec2.asuInt("PROCSEQN");
1737 int laston2 = rec2.asuInt("LASTON");
1738 if (scan1seqn == 1 && scan2seqn == 2) {
1739 if (laston1 == laston2) {
1740 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1741 //cerr<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
1742 os<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST;
1743 i +=1;
1744 }
1745 else {
1746 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1747 //cerr<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
1748 os<<LogIO::WARN<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST;
1749 }
1750 }
1751 else if (scan1seqn==2 && scan2seqn == 1) {
1752 if (laston1 == laston2) {
1753 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1754 //cerr<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<endl;
1755 os<<LogIO::WARN<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<LogIO::POST;
1756 ret = 1;
1757 break;
1758 }
1759 }
1760 else {
1761 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1762 //cerr<<"The other scan for "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<endl;
1763 os<<LogIO::WARN<<"The other scan for "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<LogIO::POST;
1764 ret = 1;
1765 break;
1766 }
1767 }
1768 }
1769 else {
1770 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1771 //cerr<<"The scan does not appear to be standard obsevation."<<endl;
1772 os<<LogIO::WARN<<"The scan does not appear to be standard obsevation."<<LogIO::POST;
1773 }
1774 //if ( i >= nscan ) break;
1775 }
1776 }
1777 else {
1778 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1779 //cerr<<"No reference to GBT_GO table."<<endl;
1780 os<<LogIO::WARN<<"No reference to GBT_GO table."<<LogIO::POST;
1781 ret = 1;
1782 }
1783 return ret;
1784}
1785
1786std::vector<double> Scantable::getDirectionVector(int whichrow) const
1787{
1788 Vector<Double> Dir = dirCol_(whichrow).getAngle("rad").getValue();
1789 std::vector<double> dir;
1790 Dir.tovector(dir);
1791 return dir;
1792}
1793
1794void asap::Scantable::reshapeSpectrum( int nmin, int nmax )
1795 throw( casa::AipsError )
1796{
1797 // assumed that all rows have same nChan
1798 Vector<Float> arr = specCol_( 0 ) ;
1799 int nChan = arr.nelements() ;
1800
1801 // if nmin < 0 or nmax < 0, nothing to do
1802 if ( nmin < 0 ) {
1803 throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
1804 }
1805 if ( nmax < 0 ) {
1806 throw( casa::indexError<int>( nmax, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
1807 }
1808
1809 // if nmin > nmax, exchange values
1810 if ( nmin > nmax ) {
1811 int tmp = nmax ;
1812 nmax = nmin ;
1813 nmin = tmp ;
1814 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1815 os << "Swap values. Applied range is ["
1816 << nmin << ", " << nmax << "]" << LogIO::POST ;
1817 }
1818
1819 // if nmin exceeds nChan, nothing to do
1820 if ( nmin >= nChan ) {
1821 throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Specified minimum exceeds nChan." ) ) ;
1822 }
1823
1824 // if nmax exceeds nChan, reset nmax to nChan
1825 if ( nmax >= nChan ) {
1826 if ( nmin == 0 ) {
1827 // nothing to do
1828 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1829 os << "Whole range is selected. Nothing to do." << LogIO::POST ;
1830 return ;
1831 }
1832 else {
1833 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1834 os << "Specified maximum exceeds nChan. Applied range is ["
1835 << nmin << ", " << nChan-1 << "]." << LogIO::POST ;
1836 nmax = nChan - 1 ;
1837 }
1838 }
1839
1840 // reshape specCol_ and flagCol_
1841 for ( int irow = 0 ; irow < nrow() ; irow++ ) {
1842 reshapeSpectrum( nmin, nmax, irow ) ;
1843 }
1844
1845 // update FREQUENCIES subtable
1846 Double refpix ;
1847 Double refval ;
1848 Double increment ;
1849 int freqnrow = freqTable_.table().nrow() ;
1850 Vector<uInt> oldId( freqnrow ) ;
1851 Vector<uInt> newId( freqnrow ) ;
1852 for ( int irow = 0 ; irow < freqnrow ; irow++ ) {
1853 freqTable_.getEntry( refpix, refval, increment, irow ) ;
1854 /***
1855 * need to shift refpix to nmin
1856 * note that channel nmin in old index will be channel 0 in new one
1857 ***/
1858 refval = refval - ( refpix - nmin ) * increment ;
1859 refpix = 0 ;
1860 freqTable_.setEntry( refpix, refval, increment, irow ) ;
1861 }
1862
1863 // update nchan
1864 int newsize = nmax - nmin + 1 ;
1865 table_.rwKeywordSet().define( "nChan", newsize ) ;
1866
1867 // update bandwidth
1868 // assumed all spectra in the scantable have same bandwidth
1869 table_.rwKeywordSet().define( "Bandwidth", increment * newsize ) ;
1870
1871 return ;
1872}
1873
1874void asap::Scantable::reshapeSpectrum( int nmin, int nmax, int irow )
1875{
1876 // reshape specCol_ and flagCol_
1877 Vector<Float> oldspec = specCol_( irow ) ;
1878 Vector<uChar> oldflag = flagsCol_( irow ) ;
1879 uInt newsize = nmax - nmin + 1 ;
1880 specCol_.put( irow, oldspec( Slice( nmin, newsize, 1 ) ) ) ;
1881 flagsCol_.put( irow, oldflag( Slice( nmin, newsize, 1 ) ) ) ;
1882
1883 return ;
1884}
1885
1886void asap::Scantable::regridChannel( int nChan, double dnu )
1887{
1888 LogIO os( LogOrigin( "Scantable", "regridChannel()", WHERE ) ) ;
1889 os << "Regrid abcissa with channel number " << nChan << " and spectral resoultion " << dnu << "Hz." << LogIO::POST ;
1890 // assumed that all rows have same nChan
1891 Vector<Float> arr = specCol_( 0 ) ;
1892 int oldsize = arr.nelements() ;
1893
1894 // if oldsize == nChan, nothing to do
1895 if ( oldsize == nChan ) {
1896 os << "Specified channel number is same as current one. Nothing to do." << LogIO::POST ;
1897 return ;
1898 }
1899
1900 // if oldChan < nChan, unphysical operation
1901 if ( oldsize < nChan ) {
1902 os << "Unphysical operation. Nothing to do." << LogIO::POST ;
1903 return ;
1904 }
1905
1906 // change channel number for specCol_ and flagCol_
1907 Vector<Float> newspec( nChan, 0 ) ;
1908 Vector<uChar> newflag( nChan, false ) ;
1909 vector<string> coordinfo = getCoordInfo() ;
1910 string oldinfo = coordinfo[0] ;
1911 coordinfo[0] = "Hz" ;
1912 setCoordInfo( coordinfo ) ;
1913 for ( int irow = 0 ; irow < nrow() ; irow++ ) {
1914 regridChannel( nChan, dnu, irow ) ;
1915 }
1916 coordinfo[0] = oldinfo ;
1917 setCoordInfo( coordinfo ) ;
1918
1919
1920 // NOTE: this method does not update metadata such as
1921 // FREQUENCIES subtable, nChan, Bandwidth, etc.
1922
1923 return ;
1924}
1925
1926void asap::Scantable::regridChannel( int nChan, double dnu, int irow )
1927{
1928 // logging
1929 //ofstream ofs( "average.log", std::ios::out | std::ios::app ) ;
1930 //ofs << "IFNO = " << getIF( irow ) << " irow = " << irow << endl ;
1931
1932 Vector<Float> oldspec = specCol_( irow ) ;
1933 Vector<uChar> oldflag = flagsCol_( irow ) ;
1934 Vector<Float> newspec( nChan, 0 ) ;
1935 Vector<uChar> newflag( nChan, false ) ;
1936
1937 // regrid
1938 vector<double> abcissa = getAbcissa( irow ) ;
1939 int oldsize = abcissa.size() ;
1940 double olddnu = abcissa[1] - abcissa[0] ;
1941 //int refChan = 0 ;
1942 //double frac = 0.0 ;
1943 //double wedge = 0.0 ;
1944 //double pile = 0.0 ;
1945 int ichan = 0 ;
1946 double wsum = 0.0 ;
1947 Vector<Float> zi( nChan+1 ) ;
1948 Vector<Float> yi( oldsize + 1 ) ;
1949 zi[0] = abcissa[0] - 0.5 * olddnu ;
1950 zi[1] = zi[1] + dnu ;
1951 for ( int ii = 2 ; ii < nChan ; ii++ )
1952 zi[ii] = zi[0] + dnu * ii ;
1953 zi[nChan] = zi[nChan-1] + dnu ;
1954 yi[0] = abcissa[0] - 0.5 * olddnu ;
1955 yi[1] = abcissa[1] + 0.5 * olddnu ;
1956 for ( int ii = 2 ; ii < oldsize ; ii++ )
1957 yi[ii] = abcissa[ii-1] + olddnu ;
1958 yi[oldsize] = abcissa[oldsize-1] + 0.5 * olddnu ;
1959 if ( dnu > 0.0 ) {
1960 for ( int ii = 0 ; ii < nChan ; ii++ ) {
1961 double zl = zi[ii] ;
1962 double zr = zi[ii+1] ;
1963 for ( int j = ichan ; j < oldsize ; j++ ) {
1964 double yl = yi[j] ;
1965 double yr = yi[j+1] ;
1966 if ( yl <= zl ) {
1967 if ( yr <= zl ) {
1968 continue ;
1969 }
1970 else if ( yr <= zr ) {
1971 newspec[ii] += oldspec[j] * ( yr - zl ) ;
1972 newflag[ii] = newflag[ii] || oldflag[j] ;
1973 wsum += ( yr - zl ) ;
1974 }
1975 else {
1976 newspec[ii] += oldspec[j] * dnu ;
1977 newflag[ii] = newflag[ii] || oldflag[j] ;
1978 wsum += dnu ;
1979 ichan = j ;
1980 break ;
1981 }
1982 }
1983 else if ( yl < zr ) {
1984 if ( yr <= zr ) {
1985 newspec[ii] += oldspec[j] * ( yr - yl ) ;
1986 newflag[ii] = newflag[ii] || oldflag[j] ;
1987 wsum += ( yr - yl ) ;
1988 }
1989 else {
1990 newspec[ii] += oldspec[j] * ( zr - yl ) ;
1991 newflag[ii] = newflag[ii] || oldflag[j] ;
1992 wsum += ( zr - yl ) ;
1993 ichan = j ;
1994 break ;
1995 }
1996 }
1997 else {
1998 ichan = j - 1 ;
1999 break ;
2000 }
2001 }
2002 if ( wsum != 0.0 )
2003 newspec[ii] /= wsum ;
2004 wsum = 0.0 ;
2005 }
2006 }
2007 else if ( dnu < 0.0 ) {
2008 for ( int ii = 0 ; ii < nChan ; ii++ ) {
2009 double zl = zi[ii] ;
2010 double zr = zi[ii+1] ;
2011 for ( int j = ichan ; j < oldsize ; j++ ) {
2012 double yl = yi[j] ;
2013 double yr = yi[j+1] ;
2014 if ( yl >= zl ) {
2015 if ( yr >= zl ) {
2016 continue ;
2017 }
2018 else if ( yr >= zr ) {
2019 newspec[ii] += oldspec[j] * abs( yr - zl ) ;
2020 newflag[ii] = newflag[ii] || oldflag[j] ;
2021 wsum += abs( yr - zl ) ;
2022 }
2023 else {
2024 newspec[ii] += oldspec[j] * abs( dnu ) ;
2025 newflag[ii] = newflag[ii] || oldflag[j] ;
2026 wsum += abs( dnu ) ;
2027 ichan = j ;
2028 break ;
2029 }
2030 }
2031 else if ( yl > zr ) {
2032 if ( yr >= zr ) {
2033 newspec[ii] += oldspec[j] * abs( yr - yl ) ;
2034 newflag[ii] = newflag[ii] || oldflag[j] ;
2035 wsum += abs( yr - yl ) ;
2036 }
2037 else {
2038 newspec[ii] += oldspec[j] * abs( zr - yl ) ;
2039 newflag[ii] = newflag[ii] || oldflag[j] ;
2040 wsum += abs( zr - yl ) ;
2041 ichan = j ;
2042 break ;
2043 }
2044 }
2045 else {
2046 ichan = j - 1 ;
2047 break ;
2048 }
2049 }
2050 if ( wsum != 0.0 )
2051 newspec[ii] /= wsum ;
2052 wsum = 0.0 ;
2053 }
2054 }
2055// * ichan = 0
2056// ***/
2057// //ofs << "olddnu = " << olddnu << ", dnu = " << dnu << endl ;
2058// pile += dnu ;
2059// wedge = olddnu * ( refChan + 1 ) ;
2060// while ( wedge < pile ) {
2061// newspec[0] += olddnu * oldspec[refChan] ;
2062// newflag[0] = newflag[0] || oldflag[refChan] ;
2063// //ofs << "channel " << refChan << " is included in new channel 0" << endl ;
2064// refChan++ ;
2065// wedge += olddnu ;
2066// wsum += olddnu ;
2067// //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
2068// }
2069// frac = ( wedge - pile ) / olddnu ;
2070// wsum += ( 1.0 - frac ) * olddnu ;
2071// newspec[0] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
2072// newflag[0] = newflag[0] || oldflag[refChan] ;
2073// //ofs << "channel " << refChan << " is partly included in new channel 0" << " with fraction of " << ( 1.0 - frac ) << endl ;
2074// //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
2075// newspec[0] /= wsum ;
2076// //ofs << "newspec[0] = " << newspec[0] << endl ;
2077// //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
2078
2079// /***
2080// * ichan = 1 - nChan-2
2081// ***/
2082// for ( int ichan = 1 ; ichan < nChan - 1 ; ichan++ ) {
2083// pile += dnu ;
2084// newspec[ichan] += frac * olddnu * oldspec[refChan] ;
2085// newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
2086// //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << frac << endl ;
2087// refChan++ ;
2088// wedge += olddnu ;
2089// wsum = frac * olddnu ;
2090// //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
2091// while ( wedge < pile ) {
2092// newspec[ichan] += olddnu * oldspec[refChan] ;
2093// newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
2094// //ofs << "channel " << refChan << " is included in new channel " << ichan << endl ;
2095// refChan++ ;
2096// wedge += olddnu ;
2097// wsum += olddnu ;
2098// //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
2099// }
2100// frac = ( wedge - pile ) / olddnu ;
2101// wsum += ( 1.0 - frac ) * olddnu ;
2102// newspec[ichan] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
2103// newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
2104// //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << ( 1.0 - frac ) << endl ;
2105// //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
2106// //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
2107// newspec[ichan] /= wsum ;
2108// //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << endl ;
2109// }
2110
2111// /***
2112// * ichan = nChan-1
2113// ***/
2114// // NOTE: Assumed that all spectra have the same bandwidth
2115// pile += dnu ;
2116// newspec[nChan-1] += frac * olddnu * oldspec[refChan] ;
2117// newflag[nChan-1] = newflag[nChan-1] || oldflag[refChan] ;
2118// //ofs << "channel " << refChan << " is partly included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
2119// refChan++ ;
2120// wedge += olddnu ;
2121// wsum = frac * olddnu ;
2122// //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
2123// for ( int jchan = refChan ; jchan < oldsize ; jchan++ ) {
2124// newspec[nChan-1] += olddnu * oldspec[jchan] ;
2125// newflag[nChan-1] = newflag[nChan-1] || oldflag[jchan] ;
2126// wsum += olddnu ;
2127// //ofs << "channel " << jchan << " is included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
2128// //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
2129// }
2130// //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
2131// //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
2132// newspec[nChan-1] /= wsum ;
2133// //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << endl ;
2134
2135// // ofs.close() ;
2136
2137 specCol_.put( irow, newspec ) ;
2138 flagsCol_.put( irow, newflag ) ;
2139
2140 return ;
2141}
2142
2143std::vector<float> Scantable::getWeather(int whichrow) const
2144{
2145 std::vector<float> out(5);
2146 //Float temperature, pressure, humidity, windspeed, windaz;
2147 weatherTable_.getEntry(out[0], out[1], out[2], out[3], out[4],
2148 mweatheridCol_(uInt(whichrow)));
2149
2150
2151 return out;
2152}
2153
2154bool Scantable::getFlagtraFast(uInt whichrow)
2155{
2156 uChar flag;
2157 Vector<uChar> flags;
2158 flagsCol_.get(whichrow, flags);
2159 flag = flags[0];
2160 for (uInt i = 1; i < flags.size(); ++i) {
2161 flag &= flags[i];
2162 }
2163 return ((flag >> 7) == 1);
2164}
2165
2166void Scantable::polyBaseline(const std::vector<bool>& mask, int order, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
2167{
2168 try {
2169 ofstream ofs;
2170 String coordInfo = "";
2171 bool hasSameNchan = true;
2172 bool outTextFile = false;
2173
2174 if (blfile != "") {
2175 ofs.open(blfile.c_str(), ios::out | ios::app);
2176 if (ofs) outTextFile = true;
2177 }
2178
2179 if (outLogger || outTextFile) {
2180 coordInfo = getCoordInfo()[0];
2181 if (coordInfo == "") coordInfo = "channel";
2182 hasSameNchan = hasSameNchanOverIFs();
2183 }
2184
2185 Fitter fitter = Fitter();
2186 fitter.setExpression("poly", order);
2187 //fitter.setIterClipping(thresClip, nIterClip);
2188
2189 int nRow = nrow();
2190 std::vector<bool> chanMask;
2191 bool showProgress;
2192 int minNRow;
2193 parseProgressInfo(progressInfo, showProgress, minNRow);
2194
2195 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
2196 chanMask = getCompositeChanMask(whichrow, mask);
2197 fitBaseline(chanMask, whichrow, fitter);
2198 setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
2199 outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "polyBaseline()", fitter);
2200 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
2201 }
2202
2203 if (outTextFile) ofs.close();
2204
2205 } catch (...) {
2206 throw;
2207 }
2208}
2209
2210void Scantable::autoPolyBaseline(const std::vector<bool>& mask, int order, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
2211{
2212 try {
2213 ofstream ofs;
2214 String coordInfo = "";
2215 bool hasSameNchan = true;
2216 bool outTextFile = false;
2217
2218 if (blfile != "") {
2219 ofs.open(blfile.c_str(), ios::out | ios::app);
2220 if (ofs) outTextFile = true;
2221 }
2222
2223 if (outLogger || outTextFile) {
2224 coordInfo = getCoordInfo()[0];
2225 if (coordInfo == "") coordInfo = "channel";
2226 hasSameNchan = hasSameNchanOverIFs();
2227 }
2228
2229 Fitter fitter = Fitter();
2230 fitter.setExpression("poly", order);
2231 //fitter.setIterClipping(thresClip, nIterClip);
2232
2233 int nRow = nrow();
2234 std::vector<bool> chanMask;
2235 int minEdgeSize = getIFNos().size()*2;
2236 STLineFinder lineFinder = STLineFinder();
2237 lineFinder.setOptions(threshold, 3, chanAvgLimit);
2238
2239 bool showProgress;
2240 int minNRow;
2241 parseProgressInfo(progressInfo, showProgress, minNRow);
2242
2243 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
2244
2245 //-------------------------------------------------------
2246 //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
2247 //-------------------------------------------------------
2248 int edgeSize = edge.size();
2249 std::vector<int> currentEdge;
2250 if (edgeSize >= 2) {
2251 int idx = 0;
2252 if (edgeSize > 2) {
2253 if (edgeSize < minEdgeSize) {
2254 throw(AipsError("Length of edge element info is less than that of IFs"));
2255 }
2256 idx = 2 * getIF(whichrow);
2257 }
2258 currentEdge.push_back(edge[idx]);
2259 currentEdge.push_back(edge[idx+1]);
2260 } else {
2261 throw(AipsError("Wrong length of edge element"));
2262 }
2263 lineFinder.setData(getSpectrum(whichrow));
2264 lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
2265 chanMask = lineFinder.getMask();
2266 //-------------------------------------------------------
2267
2268 fitBaseline(chanMask, whichrow, fitter);
2269 setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
2270
2271 outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoPolyBaseline()", fitter);
2272 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
2273 }
2274
2275 if (outTextFile) ofs.close();
2276
2277 } catch (...) {
2278 throw;
2279 }
2280}
2281
2282void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
2283{
2284 try {
2285 ofstream ofs;
2286 String coordInfo = "";
2287 bool hasSameNchan = true;
2288 bool outTextFile = false;
2289
2290 if (blfile != "") {
2291 ofs.open(blfile.c_str(), ios::out | ios::app);
2292 if (ofs) outTextFile = true;
2293 }
2294
2295 if (outLogger || outTextFile) {
2296 coordInfo = getCoordInfo()[0];
2297 if (coordInfo == "") coordInfo = "channel";
2298 hasSameNchan = hasSameNchanOverIFs();
2299 }
2300
2301 //Fitter fitter = Fitter();
2302 //fitter.setExpression("cspline", nPiece);
2303 //fitter.setIterClipping(thresClip, nIterClip);
2304
2305 int nRow = nrow();
2306 std::vector<bool> chanMask;
2307 bool showProgress;
2308 int minNRow;
2309 parseProgressInfo(progressInfo, showProgress, minNRow);
2310
2311 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
2312 chanMask = getCompositeChanMask(whichrow, mask);
2313 //fitBaseline(chanMask, whichrow, fitter);
2314 //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
2315 std::vector<int> pieceEdges;
2316 std::vector<float> params;
2317 int nClipped = 0;
2318 std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, nClipped, thresClip, nIterClip, getResidual);
2319 setSpectrum(res, whichrow);
2320 //
2321
2322 outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceEdges, params, nClipped);
2323 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
2324 }
2325
2326 if (outTextFile) ofs.close();
2327
2328 } catch (...) {
2329 throw;
2330 }
2331}
2332
2333void Scantable::autoCubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
2334{
2335 try {
2336 ofstream ofs;
2337 String coordInfo = "";
2338 bool hasSameNchan = true;
2339 bool outTextFile = false;
2340
2341 if (blfile != "") {
2342 ofs.open(blfile.c_str(), ios::out | ios::app);
2343 if (ofs) outTextFile = true;
2344 }
2345
2346 if (outLogger || outTextFile) {
2347 coordInfo = getCoordInfo()[0];
2348 if (coordInfo == "") coordInfo = "channel";
2349 hasSameNchan = hasSameNchanOverIFs();
2350 }
2351
2352 //Fitter fitter = Fitter();
2353 //fitter.setExpression("cspline", nPiece);
2354 //fitter.setIterClipping(thresClip, nIterClip);
2355
2356 int nRow = nrow();
2357 std::vector<bool> chanMask;
2358 int minEdgeSize = getIFNos().size()*2;
2359 STLineFinder lineFinder = STLineFinder();
2360 lineFinder.setOptions(threshold, 3, chanAvgLimit);
2361
2362 bool showProgress;
2363 int minNRow;
2364 parseProgressInfo(progressInfo, showProgress, minNRow);
2365
2366 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
2367
2368 //-------------------------------------------------------
2369 //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
2370 //-------------------------------------------------------
2371 int edgeSize = edge.size();
2372 std::vector<int> currentEdge;
2373 if (edgeSize >= 2) {
2374 int idx = 0;
2375 if (edgeSize > 2) {
2376 if (edgeSize < minEdgeSize) {
2377 throw(AipsError("Length of edge element info is less than that of IFs"));
2378 }
2379 idx = 2 * getIF(whichrow);
2380 }
2381 currentEdge.push_back(edge[idx]);
2382 currentEdge.push_back(edge[idx+1]);
2383 } else {
2384 throw(AipsError("Wrong length of edge element"));
2385 }
2386 lineFinder.setData(getSpectrum(whichrow));
2387 lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
2388 chanMask = lineFinder.getMask();
2389 //-------------------------------------------------------
2390
2391
2392 //fitBaseline(chanMask, whichrow, fitter);
2393 //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
2394 std::vector<int> pieceEdges;
2395 std::vector<float> params;
2396 int nClipped = 0;
2397 std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, nClipped, thresClip, nIterClip, getResidual);
2398 setSpectrum(res, whichrow);
2399 //
2400
2401 outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceEdges, params, nClipped);
2402 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
2403 }
2404
2405 if (outTextFile) ofs.close();
2406
2407 } catch (...) {
2408 throw;
2409 }
2410}
2411
2412std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, int nPiece, std::vector<int>& idxEdge, std::vector<float>& params, int& nClipped, float thresClip, int nIterClip, bool getResidual)
2413{
2414 if (data.size() != mask.size()) {
2415 throw(AipsError("data and mask sizes are not identical"));
2416 }
2417 if (nPiece < 1) {
2418 throw(AipsError("number of the sections must be one or more"));
2419 }
2420
2421 int nChan = data.size();
2422 std::vector<int> maskArray;
2423 std::vector<int> x;
2424 for (int i = 0; i < nChan; ++i) {
2425 maskArray.push_back(mask[i] ? 1 : 0);
2426 if (mask[i]) {
2427 x.push_back(i);
2428 }
2429 }
2430
2431 int initNData = x.size();
2432 if (initNData < nPiece) {
2433 throw(AipsError("too few non-flagged channels"));
2434 }
2435
2436 int nElement = (int)(floor(floor((double)(initNData/nPiece))+0.5));
2437 std::vector<double> invEdge;
2438 idxEdge.clear();
2439 idxEdge.push_back(x[0]);
2440 for (int i = 1; i < nPiece; ++i) {
2441 int valX = x[nElement*i];
2442 idxEdge.push_back(valX);
2443 invEdge.push_back(1.0/(double)valX);
2444 }
2445 idxEdge.push_back(x[x.size()-1]+1);
2446
2447 int nData = initNData;
2448 int nDOF = nPiece + 3; //number of parameters to solve, namely, 4+(nPiece-1).
2449
2450 std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1, residual;
2451 for (int i = 0; i < nChan; ++i) {
2452 double di = (double)i;
2453 double dD = (double)data[i];
2454 x1.push_back(di);
2455 x2.push_back(di*di);
2456 x3.push_back(di*di*di);
2457 z1.push_back(dD);
2458 x1z1.push_back(dD*di);
2459 x2z1.push_back(dD*di*di);
2460 x3z1.push_back(dD*di*di*di);
2461 r1.push_back(0.0);
2462 residual.push_back(0.0);
2463 }
2464
2465 for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
2466 // xMatrix : horizontal concatenation of
2467 // the least-sq. matrix (left) and an
2468 // identity matrix (right).
2469 // the right part is used to calculate the inverse matrix of the left part.
2470 double xMatrix[nDOF][2*nDOF];
2471 double zMatrix[nDOF];
2472 for (int i = 0; i < nDOF; ++i) {
2473 for (int j = 0; j < 2*nDOF; ++j) {
2474 xMatrix[i][j] = 0.0;
2475 }
2476 xMatrix[i][nDOF+i] = 1.0;
2477 zMatrix[i] = 0.0;
2478 }
2479
2480 for (int n = 0; n < nPiece; ++n) {
2481 int nUseDataInPiece = 0;
2482 for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
2483
2484 if (maskArray[i] == 0) continue;
2485
2486 xMatrix[0][0] += 1.0;
2487 xMatrix[0][1] += x1[i];
2488 xMatrix[0][2] += x2[i];
2489 xMatrix[0][3] += x3[i];
2490 xMatrix[1][1] += x2[i];
2491 xMatrix[1][2] += x3[i];
2492 xMatrix[1][3] += x2[i]*x2[i];
2493 xMatrix[2][2] += x2[i]*x2[i];
2494 xMatrix[2][3] += x3[i]*x2[i];
2495 xMatrix[3][3] += x3[i]*x3[i];
2496 zMatrix[0] += z1[i];
2497 zMatrix[1] += x1z1[i];
2498 zMatrix[2] += x2z1[i];
2499 zMatrix[3] += x3z1[i];
2500
2501 for (int j = 0; j < n; ++j) {
2502 double q = 1.0 - x1[i]*invEdge[j];
2503 q = q*q*q;
2504 xMatrix[0][j+4] += q;
2505 xMatrix[1][j+4] += q*x1[i];
2506 xMatrix[2][j+4] += q*x2[i];
2507 xMatrix[3][j+4] += q*x3[i];
2508 for (int k = 0; k < j; ++k) {
2509 double r = 1.0 - x1[i]*invEdge[k];
2510 r = r*r*r;
2511 xMatrix[k+4][j+4] += r*q;
2512 }
2513 xMatrix[j+4][j+4] += q*q;
2514 zMatrix[j+4] += q*z1[i];
2515 }
2516
2517 nUseDataInPiece++;
2518 }
2519
2520 if (nUseDataInPiece < 1) {
2521 std::vector<string> suffixOfPieceNumber(4);
2522 suffixOfPieceNumber[0] = "th";
2523 suffixOfPieceNumber[1] = "st";
2524 suffixOfPieceNumber[2] = "nd";
2525 suffixOfPieceNumber[3] = "rd";
2526 int idxNoDataPiece = (n % 10 <= 3) ? n : 0;
2527 ostringstream oss;
2528 oss << "all channels clipped or masked in " << n << suffixOfPieceNumber[idxNoDataPiece];
2529 oss << " piece of the spectrum. can't execute fitting anymore.";
2530 throw(AipsError(String(oss)));
2531 }
2532 }
2533
2534 for (int i = 0; i < nDOF; ++i) {
2535 for (int j = 0; j < i; ++j) {
2536 xMatrix[i][j] = xMatrix[j][i];
2537 }
2538 }
2539
2540 std::vector<double> invDiag;
2541 for (int i = 0; i < nDOF; ++i) {
2542 invDiag.push_back(1.0/xMatrix[i][i]);
2543 for (int j = 0; j < nDOF; ++j) {
2544 xMatrix[i][j] *= invDiag[i];
2545 }
2546 }
2547
2548 for (int k = 0; k < nDOF; ++k) {
2549 for (int i = 0; i < nDOF; ++i) {
2550 if (i != k) {
2551 double factor1 = xMatrix[k][k];
2552 double factor2 = xMatrix[i][k];
2553 for (int j = k; j < 2*nDOF; ++j) {
2554 xMatrix[i][j] *= factor1;
2555 xMatrix[i][j] -= xMatrix[k][j]*factor2;
2556 xMatrix[i][j] /= factor1;
2557 }
2558 }
2559 }
2560 double xDiag = xMatrix[k][k];
2561 for (int j = k; j < 2*nDOF; ++j) {
2562 xMatrix[k][j] /= xDiag;
2563 }
2564 }
2565
2566 for (int i = 0; i < nDOF; ++i) {
2567 for (int j = 0; j < nDOF; ++j) {
2568 xMatrix[i][nDOF+j] *= invDiag[j];
2569 }
2570 }
2571 //compute a vector y which consists of the coefficients of the best-fit spline curves
2572 //(a0,a1,a2,a3(,b3,c3,...)), namely, the ones for the leftmost piece and the ones of
2573 //cubic terms for the other pieces (in case nPiece>1).
2574 std::vector<double> y;
2575 y.clear();
2576 for (int i = 0; i < nDOF; ++i) {
2577 y.push_back(0.0);
2578 for (int j = 0; j < nDOF; ++j) {
2579 y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
2580 }
2581 }
2582
2583 double a0 = y[0];
2584 double a1 = y[1];
2585 double a2 = y[2];
2586 double a3 = y[3];
2587 params.clear();
2588
2589 for (int n = 0; n < nPiece; ++n) {
2590 for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
2591 r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
2592 residual[i] = z1[i] - r1[i];
2593 }
2594 params.push_back(a0);
2595 params.push_back(a1);
2596 params.push_back(a2);
2597 params.push_back(a3);
2598
2599 if (n == nPiece-1) break;
2600
2601 double d = y[4+n];
2602 double iE = invEdge[n];
2603 a0 += d;
2604 a1 -= 3.0*d*iE;
2605 a2 += 3.0*d*iE*iE;
2606 a3 -= d*iE*iE*iE;
2607 }
2608
2609 if ((nClip == nIterClip) || (thresClip <= 0.0)) {
2610 break;
2611 } else {
2612 double stdDev = 0.0;
2613 for (int i = 0; i < nChan; ++i) {
2614 stdDev += residual[i]*residual[i]*(double)maskArray[i];
2615 }
2616 stdDev = sqrt(stdDev/(double)nData);
2617
2618 double thres = stdDev * thresClip;
2619 int newNData = 0;
2620 for (int i = 0; i < nChan; ++i) {
2621 if (abs(residual[i]) >= thres) {
2622 maskArray[i] = 0;
2623 }
2624 if (maskArray[i] > 0) {
2625 newNData++;
2626 }
2627 }
2628 if (newNData == nData) {
2629 break; //no more flag to add. iteration stops.
2630 } else {
2631 nData = newNData;
2632 }
2633 }
2634 }
2635
2636 nClipped = initNData - nData;
2637
2638 std::vector<float> result;
2639 if (getResidual) {
2640 for (int i = 0; i < nChan; ++i) {
2641 result.push_back((float)residual[i]);
2642 }
2643 } else {
2644 for (int i = 0; i < nChan; ++i) {
2645 result.push_back((float)r1[i]);
2646 }
2647 }
2648
2649 return result;
2650}
2651
2652 void Scantable::selectWaveNumbers(const int whichrow, const std::vector<bool>& chanMask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
2653{
2654 nWaves.clear();
2655
2656 if (applyFFT) {
2657 string fftThAttr;
2658 float fftThSigma;
2659 int fftThTop;
2660 parseThresholdExpression(fftThresh, fftThAttr, fftThSigma, fftThTop);
2661 doSelectWaveNumbers(whichrow, chanMask, fftMethod, fftThSigma, fftThTop, fftThAttr, nWaves);
2662 }
2663
2664 addAuxWaveNumbers(addNWaves, rejectNWaves, nWaves);
2665}
2666
2667void Scantable::parseThresholdExpression(const std::string& fftThresh, std::string& fftThAttr, float& fftThSigma, int& fftThTop)
2668{
2669 uInt idxSigma = fftThresh.find("sigma");
2670 uInt idxTop = fftThresh.find("top");
2671
2672 if (idxSigma == fftThresh.size() - 5) {
2673 std::istringstream is(fftThresh.substr(0, fftThresh.size() - 5));
2674 is >> fftThSigma;
2675 fftThAttr = "sigma";
2676 } else if (idxTop == 0) {
2677 std::istringstream is(fftThresh.substr(3));
2678 is >> fftThTop;
2679 fftThAttr = "top";
2680 } else {
2681 bool isNumber = true;
2682 for (uInt i = 0; i < fftThresh.size()-1; ++i) {
2683 char ch = (fftThresh.substr(i, 1).c_str())[0];
2684 if (!(isdigit(ch) || (fftThresh.substr(i, 1) == "."))) {
2685 isNumber = false;
2686 break;
2687 }
2688 }
2689 if (isNumber) {
2690 std::istringstream is(fftThresh);
2691 is >> fftThSigma;
2692 fftThAttr = "sigma";
2693 } else {
2694 throw(AipsError("fftthresh has a wrong value"));
2695 }
2696 }
2697}
2698
2699void 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)
2700{
2701 std::vector<float> fspec;
2702 if (fftMethod == "fft") {
2703 fspec = execFFT(whichrow, chanMask, false, true);
2704 //} else if (fftMethod == "lsp") {
2705 // fspec = lombScarglePeriodogram(whichrow);
2706 }
2707
2708 if (fftThAttr == "sigma") {
2709 float mean = 0.0;
2710 float mean2 = 0.0;
2711 for (uInt i = 0; i < fspec.size(); ++i) {
2712 mean += fspec[i];
2713 mean2 += fspec[i]*fspec[i];
2714 }
2715 mean /= float(fspec.size());
2716 mean2 /= float(fspec.size());
2717 float thres = mean + fftThSigma * float(sqrt(mean2 - mean*mean));
2718
2719 for (uInt i = 0; i < fspec.size(); ++i) {
2720 if (fspec[i] >= thres) {
2721 nWaves.push_back(i);
2722 }
2723 }
2724
2725 } else if (fftThAttr == "top") {
2726 for (int i = 0; i < fftThTop; ++i) {
2727 float max = 0.0;
2728 int maxIdx = 0;
2729 for (uInt j = 0; j < fspec.size(); ++j) {
2730 if (fspec[j] > max) {
2731 max = fspec[j];
2732 maxIdx = j;
2733 }
2734 }
2735 nWaves.push_back(maxIdx);
2736 fspec[maxIdx] = 0.0;
2737 }
2738
2739 }
2740
2741 if (nWaves.size() > 1) {
2742 sort(nWaves.begin(), nWaves.end());
2743 }
2744}
2745
2746void Scantable::addAuxWaveNumbers(const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
2747{
2748 for (uInt i = 0; i < addNWaves.size(); ++i) {
2749 bool found = false;
2750 for (uInt j = 0; j < nWaves.size(); ++j) {
2751 if (nWaves[j] == addNWaves[i]) {
2752 found = true;
2753 break;
2754 }
2755 }
2756 if (!found) nWaves.push_back(addNWaves[i]);
2757 }
2758
2759 for (uInt i = 0; i < rejectNWaves.size(); ++i) {
2760 for (std::vector<int>::iterator j = nWaves.begin(); j != nWaves.end(); ) {
2761 if (*j == rejectNWaves[i]) {
2762 j = nWaves.erase(j);
2763 } else {
2764 ++j;
2765 }
2766 }
2767 }
2768
2769 if (nWaves.size() > 1) {
2770 sort(nWaves.begin(), nWaves.end());
2771 unique(nWaves.begin(), nWaves.end());
2772 }
2773}
2774
2775void Scantable::sinusoidBaseline(const std::vector<bool>& mask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, float thresClip, int nIterClip, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
2776{
2777 try {
2778 ofstream ofs;
2779 String coordInfo = "";
2780 bool hasSameNchan = true;
2781 bool outTextFile = false;
2782
2783 if (blfile != "") {
2784 ofs.open(blfile.c_str(), ios::out | ios::app);
2785 if (ofs) outTextFile = true;
2786 }
2787
2788 if (outLogger || outTextFile) {
2789 coordInfo = getCoordInfo()[0];
2790 if (coordInfo == "") coordInfo = "channel";
2791 hasSameNchan = hasSameNchanOverIFs();
2792 }
2793
2794 //Fitter fitter = Fitter();
2795 //fitter.setExpression("sinusoid", nWaves);
2796 //fitter.setIterClipping(thresClip, nIterClip);
2797
2798 int nRow = nrow();
2799 std::vector<bool> chanMask;
2800 std::vector<int> nWaves;
2801
2802 bool showProgress;
2803 int minNRow;
2804 parseProgressInfo(progressInfo, showProgress, minNRow);
2805
2806 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
2807 chanMask = getCompositeChanMask(whichrow, mask);
2808 selectWaveNumbers(whichrow, chanMask, applyFFT, fftMethod, fftThresh, addNWaves, rejectNWaves, nWaves);
2809
2810 //FOR DEBUGGING------------
2811 if (whichrow < 0) {// == nRow -1) {
2812 cout << "+++ i=" << setw(3) << whichrow << ", IF=" << setw(2) << getIF(whichrow);
2813 if (applyFFT) {
2814 cout << "[ ";
2815 for (uInt j = 0; j < nWaves.size(); ++j) {
2816 cout << nWaves[j] << ", ";
2817 }
2818 cout << " ] " << endl;
2819 }
2820 cout << flush;
2821 }
2822 //-------------------------
2823
2824 //fitBaseline(chanMask, whichrow, fitter);
2825 //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
2826 std::vector<float> params;
2827 int nClipped = 0;
2828 std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, nClipped, thresClip, nIterClip, getResidual);
2829 setSpectrum(res, whichrow);
2830 //
2831
2832 outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", params, nClipped);
2833 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
2834 }
2835
2836 if (outTextFile) ofs.close();
2837
2838 } catch (...) {
2839 throw;
2840 }
2841}
2842
2843void Scantable::autoSinusoidBaseline(const std::vector<bool>& mask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
2844{
2845 try {
2846 ofstream ofs;
2847 String coordInfo = "";
2848 bool hasSameNchan = true;
2849 bool outTextFile = false;
2850
2851 if (blfile != "") {
2852 ofs.open(blfile.c_str(), ios::out | ios::app);
2853 if (ofs) outTextFile = true;
2854 }
2855
2856 if (outLogger || outTextFile) {
2857 coordInfo = getCoordInfo()[0];
2858 if (coordInfo == "") coordInfo = "channel";
2859 hasSameNchan = hasSameNchanOverIFs();
2860 }
2861
2862 //Fitter fitter = Fitter();
2863 //fitter.setExpression("sinusoid", nWaves);
2864 //fitter.setIterClipping(thresClip, nIterClip);
2865
2866 int nRow = nrow();
2867 std::vector<bool> chanMask;
2868 std::vector<int> nWaves;
2869
2870 int minEdgeSize = getIFNos().size()*2;
2871 STLineFinder lineFinder = STLineFinder();
2872 lineFinder.setOptions(threshold, 3, chanAvgLimit);
2873
2874 bool showProgress;
2875 int minNRow;
2876 parseProgressInfo(progressInfo, showProgress, minNRow);
2877
2878 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
2879
2880 //-------------------------------------------------------
2881 //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
2882 //-------------------------------------------------------
2883 int edgeSize = edge.size();
2884 std::vector<int> currentEdge;
2885 if (edgeSize >= 2) {
2886 int idx = 0;
2887 if (edgeSize > 2) {
2888 if (edgeSize < minEdgeSize) {
2889 throw(AipsError("Length of edge element info is less than that of IFs"));
2890 }
2891 idx = 2 * getIF(whichrow);
2892 }
2893 currentEdge.push_back(edge[idx]);
2894 currentEdge.push_back(edge[idx+1]);
2895 } else {
2896 throw(AipsError("Wrong length of edge element"));
2897 }
2898 lineFinder.setData(getSpectrum(whichrow));
2899 lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
2900 chanMask = lineFinder.getMask();
2901 //-------------------------------------------------------
2902
2903 selectWaveNumbers(whichrow, chanMask, applyFFT, fftMethod, fftThresh, addNWaves, rejectNWaves, nWaves);
2904
2905 //fitBaseline(chanMask, whichrow, fitter);
2906 //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
2907 std::vector<float> params;
2908 int nClipped = 0;
2909 std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, nClipped, thresClip, nIterClip, getResidual);
2910 setSpectrum(res, whichrow);
2911 //
2912
2913 outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", params, nClipped);
2914 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
2915 }
2916
2917 if (outTextFile) ofs.close();
2918
2919 } catch (...) {
2920 throw;
2921 }
2922}
2923
2924std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data, const std::vector<bool>& mask, const std::vector<int>& waveNumbers, std::vector<float>& params, int& nClipped, float thresClip, int nIterClip, bool getResidual)
2925{
2926 if (data.size() != mask.size()) {
2927 throw(AipsError("data and mask sizes are not identical"));
2928 }
2929 if (data.size() < 2) {
2930 throw(AipsError("data size is too short"));
2931 }
2932 if (waveNumbers.size() == 0) {
2933 throw(AipsError("no wave numbers given"));
2934 }
2935 std::vector<int> nWaves; // sorted and uniqued array of wave numbers
2936 nWaves.reserve(waveNumbers.size());
2937 copy(waveNumbers.begin(), waveNumbers.end(), back_inserter(nWaves));
2938 sort(nWaves.begin(), nWaves.end());
2939 std::vector<int>::iterator end_it = unique(nWaves.begin(), nWaves.end());
2940 nWaves.erase(end_it, nWaves.end());
2941
2942 int minNWaves = nWaves[0];
2943 if (minNWaves < 0) {
2944 throw(AipsError("wave number must be positive or zero (i.e. constant)"));
2945 }
2946 bool hasConstantTerm = (minNWaves == 0);
2947
2948 int nChan = data.size();
2949 std::vector<int> maskArray;
2950 std::vector<int> x;
2951 for (int i = 0; i < nChan; ++i) {
2952 maskArray.push_back(mask[i] ? 1 : 0);
2953 if (mask[i]) {
2954 x.push_back(i);
2955 }
2956 }
2957
2958 int initNData = x.size();
2959
2960 int nData = initNData;
2961 int nDOF = nWaves.size() * 2 - (hasConstantTerm ? 1 : 0); //number of parameters to solve.
2962
2963 const double PI = 6.0 * asin(0.5); // PI (= 3.141592653...)
2964 double baseXFactor = 2.0*PI/(double)(nChan-1); //the denominator (nChan-1) should be changed to (xdata[nChan-1]-xdata[0]) for accepting x-values given in velocity or frequency when this function is moved to fitter. (2011/03/30 WK)
2965
2966 // xArray : contains elemental values for computing the least-square matrix.
2967 // xArray.size() is nDOF and xArray[*].size() is nChan.
2968 // Each xArray element are as follows:
2969 // xArray[0] = {1.0, 1.0, 1.0, ..., 1.0},
2970 // xArray[2n-1] = {sin(nPI/L*x[0]), sin(nPI/L*x[1]), ..., sin(nPI/L*x[nChan])},
2971 // xArray[2n] = {cos(nPI/L*x[0]), cos(nPI/L*x[1]), ..., cos(nPI/L*x[nChan])},
2972 // where (1 <= n <= nMaxWavesInSW),
2973 // or,
2974 // xArray[2n-1] = {sin(wn[n]PI/L*x[0]), sin(wn[n]PI/L*x[1]), ..., sin(wn[n]PI/L*x[nChan])},
2975 // xArray[2n] = {cos(wn[n]PI/L*x[0]), cos(wn[n]PI/L*x[1]), ..., cos(wn[n]PI/L*x[nChan])},
2976 // where wn[n] denotes waveNumbers[n] (1 <= n <= waveNumbers.size()).
2977 std::vector<std::vector<double> > xArray;
2978 if (hasConstantTerm) {
2979 std::vector<double> xu;
2980 for (int j = 0; j < nChan; ++j) {
2981 xu.push_back(1.0);
2982 }
2983 xArray.push_back(xu);
2984 }
2985 for (uInt i = (hasConstantTerm ? 1 : 0); i < nWaves.size(); ++i) {
2986 double xFactor = baseXFactor*(double)nWaves[i];
2987 std::vector<double> xs, xc;
2988 xs.clear();
2989 xc.clear();
2990 for (int j = 0; j < nChan; ++j) {
2991 xs.push_back(sin(xFactor*(double)j));
2992 xc.push_back(cos(xFactor*(double)j));
2993 }
2994 xArray.push_back(xs);
2995 xArray.push_back(xc);
2996 }
2997
2998 std::vector<double> z1, r1, residual;
2999 for (int i = 0; i < nChan; ++i) {
3000 z1.push_back((double)data[i]);
3001 r1.push_back(0.0);
3002 residual.push_back(0.0);
3003 }
3004
3005 for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
3006 // xMatrix : horizontal concatenation of
3007 // the least-sq. matrix (left) and an
3008 // identity matrix (right).
3009 // the right part is used to calculate the inverse matrix of the left part.
3010 double xMatrix[nDOF][2*nDOF];
3011 double zMatrix[nDOF];
3012 for (int i = 0; i < nDOF; ++i) {
3013 for (int j = 0; j < 2*nDOF; ++j) {
3014 xMatrix[i][j] = 0.0;
3015 }
3016 xMatrix[i][nDOF+i] = 1.0;
3017 zMatrix[i] = 0.0;
3018 }
3019
3020 int nUseData = 0;
3021 for (int k = 0; k < nChan; ++k) {
3022 if (maskArray[k] == 0) continue;
3023
3024 for (int i = 0; i < nDOF; ++i) {
3025 for (int j = i; j < nDOF; ++j) {
3026 xMatrix[i][j] += xArray[i][k] * xArray[j][k];
3027 }
3028 zMatrix[i] += z1[k] * xArray[i][k];
3029 }
3030
3031 nUseData++;
3032 }
3033
3034 if (nUseData < 1) {
3035 throw(AipsError("all channels clipped or masked. can't execute fitting anymore."));
3036 }
3037
3038 for (int i = 0; i < nDOF; ++i) {
3039 for (int j = 0; j < i; ++j) {
3040 xMatrix[i][j] = xMatrix[j][i];
3041 }
3042 }
3043
3044 std::vector<double> invDiag;
3045 for (int i = 0; i < nDOF; ++i) {
3046 invDiag.push_back(1.0/xMatrix[i][i]);
3047 for (int j = 0; j < nDOF; ++j) {
3048 xMatrix[i][j] *= invDiag[i];
3049 }
3050 }
3051
3052 for (int k = 0; k < nDOF; ++k) {
3053 for (int i = 0; i < nDOF; ++i) {
3054 if (i != k) {
3055 double factor1 = xMatrix[k][k];
3056 double factor2 = xMatrix[i][k];
3057 for (int j = k; j < 2*nDOF; ++j) {
3058 xMatrix[i][j] *= factor1;
3059 xMatrix[i][j] -= xMatrix[k][j]*factor2;
3060 xMatrix[i][j] /= factor1;
3061 }
3062 }
3063 }
3064 double xDiag = xMatrix[k][k];
3065 for (int j = k; j < 2*nDOF; ++j) {
3066 xMatrix[k][j] /= xDiag;
3067 }
3068 }
3069
3070 for (int i = 0; i < nDOF; ++i) {
3071 for (int j = 0; j < nDOF; ++j) {
3072 xMatrix[i][nDOF+j] *= invDiag[j];
3073 }
3074 }
3075 //compute a vector y which consists of the coefficients of the sinusoids forming the
3076 //best-fit curves (a0,s1,c1,s2,c2,...), where a0 is constant and s* and c* are of sine
3077 //and cosine functions, respectively.
3078 std::vector<double> y;
3079 params.clear();
3080 for (int i = 0; i < nDOF; ++i) {
3081 y.push_back(0.0);
3082 for (int j = 0; j < nDOF; ++j) {
3083 y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
3084 }
3085 params.push_back(y[i]);
3086 }
3087
3088 for (int i = 0; i < nChan; ++i) {
3089 r1[i] = y[0];
3090 for (int j = 1; j < nDOF; ++j) {
3091 r1[i] += y[j]*xArray[j][i];
3092 }
3093 residual[i] = z1[i] - r1[i];
3094 }
3095
3096 if ((nClip == nIterClip) || (thresClip <= 0.0)) {
3097 break;
3098 } else {
3099 double stdDev = 0.0;
3100 for (int i = 0; i < nChan; ++i) {
3101 stdDev += residual[i]*residual[i]*(double)maskArray[i];
3102 }
3103 stdDev = sqrt(stdDev/(double)nData);
3104
3105 double thres = stdDev * thresClip;
3106 int newNData = 0;
3107 for (int i = 0; i < nChan; ++i) {
3108 if (abs(residual[i]) >= thres) {
3109 maskArray[i] = 0;
3110 }
3111 if (maskArray[i] > 0) {
3112 newNData++;
3113 }
3114 }
3115 if (newNData == nData) {
3116 break; //no more flag to add. iteration stops.
3117 } else {
3118 nData = newNData;
3119 }
3120 }
3121 }
3122
3123 nClipped = initNData - nData;
3124
3125 std::vector<float> result;
3126 if (getResidual) {
3127 for (int i = 0; i < nChan; ++i) {
3128 result.push_back((float)residual[i]);
3129 }
3130 } else {
3131 for (int i = 0; i < nChan; ++i) {
3132 result.push_back((float)r1[i]);
3133 }
3134 }
3135
3136 return result;
3137}
3138
3139void Scantable::fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter)
3140{
3141 std::vector<double> dAbcissa = getAbcissa(whichrow);
3142 std::vector<float> abcissa;
3143 for (uInt i = 0; i < dAbcissa.size(); ++i) {
3144 abcissa.push_back((float)dAbcissa[i]);
3145 }
3146 std::vector<float> spec = getSpectrum(whichrow);
3147
3148 fitter.setData(abcissa, spec, mask);
3149 fitter.lfit();
3150}
3151
3152std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask)
3153{
3154 std::vector<bool> mask = getMask(whichrow);
3155 uInt maskSize = mask.size();
3156 if (maskSize != inMask.size()) {
3157 throw(AipsError("mask sizes are not the same."));
3158 }
3159 for (uInt i = 0; i < maskSize; ++i) {
3160 mask[i] = mask[i] && inMask[i];
3161 }
3162
3163 return mask;
3164}
3165
3166/*
3167std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder)
3168{
3169 int edgeSize = edge.size();
3170 std::vector<int> currentEdge;
3171 if (edgeSize >= 2) {
3172 int idx = 0;
3173 if (edgeSize > 2) {
3174 if (edgeSize < minEdgeSize) {
3175 throw(AipsError("Length of edge element info is less than that of IFs"));
3176 }
3177 idx = 2 * getIF(whichrow);
3178 }
3179 currentEdge.push_back(edge[idx]);
3180 currentEdge.push_back(edge[idx+1]);
3181 } else {
3182 throw(AipsError("Wrong length of edge element"));
3183 }
3184
3185 lineFinder.setData(getSpectrum(whichrow));
3186 lineFinder.findLines(getCompositeChanMask(whichrow, inMask), currentEdge, whichrow);
3187
3188 return lineFinder.getMask();
3189}
3190*/
3191
3192/* for poly. the variations of outputFittingResult() should be merged into one eventually (2011/3/10 WK) */
3193void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, Fitter& fitter)
3194{
3195 if (outLogger || outTextFile) {
3196 std::vector<float> params = fitter.getParameters();
3197 std::vector<bool> fixed = fitter.getFixedParameters();
3198 float rms = getRms(chanMask, whichrow);
3199 String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
3200
3201 if (outLogger) {
3202 LogIO ols(LogOrigin("Scantable", funcName, WHERE));
3203 ols << formatBaselineParams(params, fixed, rms, -1, masklist, whichrow, false) << LogIO::POST ;
3204 }
3205 if (outTextFile) {
3206 ofs << formatBaselineParams(params, fixed, rms, -1, masklist, whichrow, true) << flush;
3207 }
3208 }
3209}
3210
3211/* for cspline. will be merged once cspline is available in fitter (2011/3/10 WK) */
3212void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<int>& edge, const std::vector<float>& params, const int nClipped)
3213{
3214 if (outLogger || outTextFile) {
3215 float rms = getRms(chanMask, whichrow);
3216 String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
3217 std::vector<bool> fixed;
3218 fixed.clear();
3219
3220 if (outLogger) {
3221 LogIO ols(LogOrigin("Scantable", funcName, WHERE));
3222 ols << formatPiecewiseBaselineParams(edge, params, fixed, rms, nClipped, masklist, whichrow, false) << LogIO::POST ;
3223 }
3224 if (outTextFile) {
3225 ofs << formatPiecewiseBaselineParams(edge, params, fixed, rms, nClipped, masklist, whichrow, true) << flush;
3226 }
3227 }
3228}
3229
3230/* for sinusoid. will be merged once sinusoid is available in fitter (2011/3/10 WK) */
3231void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<float>& params, const int nClipped)
3232{
3233 if (outLogger || outTextFile) {
3234 float rms = getRms(chanMask, whichrow);
3235 String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
3236 std::vector<bool> fixed;
3237 fixed.clear();
3238
3239 if (outLogger) {
3240 LogIO ols(LogOrigin("Scantable", funcName, WHERE));
3241 ols << formatBaselineParams(params, fixed, rms, nClipped, masklist, whichrow, false) << LogIO::POST ;
3242 }
3243 if (outTextFile) {
3244 ofs << formatBaselineParams(params, fixed, rms, nClipped, masklist, whichrow, true) << flush;
3245 }
3246 }
3247}
3248
3249void Scantable::parseProgressInfo(const std::string& progressInfo, bool& showProgress, int& minNRow)
3250{
3251 int idxDelimiter = progressInfo.find(",");
3252 if (idxDelimiter < 0) {
3253 throw(AipsError("wrong value in 'showprogress' parameter")) ;
3254 }
3255 showProgress = (progressInfo.substr(0, idxDelimiter) == "true");
3256 std::istringstream is(progressInfo.substr(idxDelimiter+1));
3257 is >> minNRow;
3258}
3259
3260void Scantable::showProgressOnTerminal(const int nProcessed, const int nTotal, const bool showProgress, const int nTotalThreshold)
3261{
3262 if (showProgress && (nTotal >= nTotalThreshold)) {
3263 int nInterval = int(floor(double(nTotal)/100.0));
3264 if (nInterval == 0) nInterval++;
3265
3266 if (nProcessed % nInterval == 0) {
3267 printf("\r"); //go to the head of line
3268 printf("\x1b[31m\x1b[1m"); //set red color, highlighted
3269 printf("[%3d%%]", (int)(100.0*(double(nProcessed+1))/(double(nTotal))) );
3270 printf("\x1b[39m\x1b[0m"); //set default attributes
3271 fflush(NULL);
3272 }
3273
3274 if (nProcessed == nTotal - 1) {
3275 printf("\r\x1b[K"); //clear
3276 fflush(NULL);
3277 }
3278
3279 }
3280}
3281
3282std::vector<float> Scantable::execFFT(const int whichrow, const std::vector<bool>& inMask, bool getRealImag, bool getAmplitudeOnly)
3283{
3284 std::vector<bool> mask = getMask(whichrow);
3285
3286 if (inMask.size() > 0) {
3287 uInt maskSize = mask.size();
3288 if (maskSize != inMask.size()) {
3289 throw(AipsError("mask sizes are not the same."));
3290 }
3291 for (uInt i = 0; i < maskSize; ++i) {
3292 mask[i] = mask[i] && inMask[i];
3293 }
3294 }
3295
3296 Vector<Float> spec = getSpectrum(whichrow);
3297 mathutil::doZeroOrderInterpolation(spec, mask);
3298
3299 FFTServer<Float,Complex> ffts;
3300 Vector<Complex> fftres;
3301 ffts.fft0(fftres, spec);
3302
3303 std::vector<float> res;
3304 float norm = float(2.0/double(spec.size()));
3305
3306 if (getRealImag) {
3307 for (uInt i = 0; i < fftres.size(); ++i) {
3308 res.push_back(real(fftres[i])*norm);
3309 res.push_back(imag(fftres[i])*norm);
3310 }
3311 } else {
3312 for (uInt i = 0; i < fftres.size(); ++i) {
3313 res.push_back(abs(fftres[i])*norm);
3314 if (!getAmplitudeOnly) res.push_back(arg(fftres[i]));
3315 }
3316 }
3317
3318 return res;
3319}
3320
3321
3322float Scantable::getRms(const std::vector<bool>& mask, int whichrow)
3323{
3324 Vector<Float> spec;
3325 specCol_.get(whichrow, spec);
3326
3327 float mean = 0.0;
3328 float smean = 0.0;
3329 int n = 0;
3330 for (uInt i = 0; i < spec.nelements(); ++i) {
3331 if (mask[i]) {
3332 mean += spec[i];
3333 smean += spec[i]*spec[i];
3334 n++;
3335 }
3336 }
3337
3338 mean /= (float)n;
3339 smean /= (float)n;
3340
3341 return sqrt(smean - mean*mean);
3342}
3343
3344
3345std::string Scantable::formatBaselineParamsHeader(int whichrow, const std::string& masklist, bool verbose) const
3346{
3347 ostringstream oss;
3348
3349 if (verbose) {
3350 oss << " Scan[" << getScan(whichrow) << "]";
3351 oss << " Beam[" << getBeam(whichrow) << "]";
3352 oss << " IF[" << getIF(whichrow) << "]";
3353 oss << " Pol[" << getPol(whichrow) << "]";
3354 oss << " Cycle[" << getCycle(whichrow) << "]: " << endl;
3355 oss << "Fitter range = " << masklist << endl;
3356 oss << "Baseline parameters" << endl;
3357 oss << flush;
3358 }
3359
3360 return String(oss);
3361}
3362
3363std::string Scantable::formatBaselineParamsFooter(float rms, int nClipped, bool verbose) const
3364{
3365 ostringstream oss;
3366
3367 if (verbose) {
3368 oss << "Results of baseline fit" << endl;
3369 oss << " rms = " << setprecision(6) << rms << endl;
3370 if (nClipped >= 0) {
3371 oss << " Number of clipped channels = " << nClipped << endl;
3372 }
3373 for (int i = 0; i < 60; ++i) {
3374 oss << "-";
3375 }
3376 oss << endl;
3377 oss << flush;
3378 }
3379
3380 return String(oss);
3381}
3382
3383std::string Scantable::formatBaselineParams(const std::vector<float>& params,
3384 const std::vector<bool>& fixed,
3385 float rms,
3386 int nClipped,
3387 const std::string& masklist,
3388 int whichrow,
3389 bool verbose,
3390 int start, int count,
3391 bool resetparamid) const
3392{
3393 int nParam = (int)(params.size());
3394
3395 if (nParam < 1) {
3396 return(" Not fitted");
3397 } else {
3398
3399 ostringstream oss;
3400 oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
3401
3402 if (start < 0) start = 0;
3403 if (count < 0) count = nParam;
3404 int end = start + count;
3405 if (end > nParam) end = nParam;
3406 int paramidoffset = (resetparamid) ? (-start) : 0;
3407
3408 for (int i = start; i < end; ++i) {
3409 if (i > start) {
3410 oss << ",";
3411 }
3412 std::string sFix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
3413 oss << " p" << (i+paramidoffset) << sFix << "= " << right << setw(13) << setprecision(6) << params[i];
3414 }
3415
3416 oss << endl;
3417 oss << formatBaselineParamsFooter(rms, nClipped, verbose);
3418
3419 return String(oss);
3420 }
3421
3422}
3423
3424 std::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) const
3425{
3426 int nOutParam = (int)(params.size());
3427 int nPiece = (int)(ranges.size()) - 1;
3428
3429 if (nOutParam < 1) {
3430 return(" Not fitted");
3431 } else if (nPiece < 0) {
3432 return formatBaselineParams(params, fixed, rms, nClipped, masklist, whichrow, verbose);
3433 } else if (nPiece < 1) {
3434 return(" Bad count of the piece edge info");
3435 } else if (nOutParam % nPiece != 0) {
3436 return(" Bad count of the output baseline parameters");
3437 } else {
3438
3439 int nParam = nOutParam / nPiece;
3440
3441 ostringstream oss;
3442 oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
3443
3444 stringstream ss;
3445 ss << ranges[nPiece] << flush;
3446 int wRange = ss.str().size() * 2 + 5;
3447
3448 for (int i = 0; i < nPiece; ++i) {
3449 ss.str("");
3450 ss << " [" << ranges[i] << "," << (ranges[i+1]-1) << "]";
3451 oss << left << setw(wRange) << ss.str();
3452 oss << formatBaselineParams(params, fixed, rms, 0, masklist, whichrow, false, i*nParam, nParam, true);
3453 }
3454
3455 oss << formatBaselineParamsFooter(rms, nClipped, verbose);
3456
3457 return String(oss);
3458 }
3459
3460}
3461
3462bool Scantable::hasSameNchanOverIFs()
3463{
3464 int nIF = nif(-1);
3465 int nCh;
3466 int totalPositiveNChan = 0;
3467 int nPositiveNChan = 0;
3468
3469 for (int i = 0; i < nIF; ++i) {
3470 nCh = nchan(i);
3471 if (nCh > 0) {
3472 totalPositiveNChan += nCh;
3473 nPositiveNChan++;
3474 }
3475 }
3476
3477 return (totalPositiveNChan == (nPositiveNChan * nchan(0)));
3478}
3479
3480std::string Scantable::getMaskRangeList(const std::vector<bool>& mask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, bool verbose)
3481{
3482 if (mask.size() < 2) {
3483 throw(AipsError("The mask elements should be > 1"));
3484 }
3485 int IF = getIF(whichrow);
3486 if (mask.size() != (uInt)nchan(IF)) {
3487 throw(AipsError("Number of channels in scantable != number of mask elements"));
3488 }
3489
3490 if (verbose) {
3491 LogIO logOs(LogOrigin("Scantable", "getMaskRangeList()", WHERE));
3492 logOs << LogIO::WARN << "The current mask window unit is " << coordInfo;
3493 if (!hasSameNchan) {
3494 logOs << endl << "This mask is only valid for IF=" << IF;
3495 }
3496 logOs << LogIO::POST;
3497 }
3498
3499 std::vector<double> abcissa = getAbcissa(whichrow);
3500 std::vector<int> edge = getMaskEdgeIndices(mask);
3501
3502 ostringstream oss;
3503 oss.setf(ios::fixed);
3504 oss << setprecision(1) << "[";
3505 for (uInt i = 0; i < edge.size(); i+=2) {
3506 if (i > 0) oss << ",";
3507 oss << "[" << (float)abcissa[edge[i]] << "," << (float)abcissa[edge[i+1]] << "]";
3508 }
3509 oss << "]" << flush;
3510
3511 return String(oss);
3512}
3513
3514std::vector<int> Scantable::getMaskEdgeIndices(const std::vector<bool>& mask)
3515{
3516 if (mask.size() < 2) {
3517 throw(AipsError("The mask elements should be > 1"));
3518 }
3519
3520 std::vector<int> out, startIndices, endIndices;
3521 int maskSize = mask.size();
3522
3523 startIndices.clear();
3524 endIndices.clear();
3525
3526 if (mask[0]) {
3527 startIndices.push_back(0);
3528 }
3529 for (int i = 1; i < maskSize; ++i) {
3530 if ((!mask[i-1]) && mask[i]) {
3531 startIndices.push_back(i);
3532 } else if (mask[i-1] && (!mask[i])) {
3533 endIndices.push_back(i-1);
3534 }
3535 }
3536 if (mask[maskSize-1]) {
3537 endIndices.push_back(maskSize-1);
3538 }
3539
3540 if (startIndices.size() != endIndices.size()) {
3541 throw(AipsError("Inconsistent Mask Size: bad data?"));
3542 }
3543 for (uInt i = 0; i < startIndices.size(); ++i) {
3544 if (startIndices[i] > endIndices[i]) {
3545 throw(AipsError("Mask start index > mask end index"));
3546 }
3547 }
3548
3549 out.clear();
3550 for (uInt i = 0; i < startIndices.size(); ++i) {
3551 out.push_back(startIndices[i]);
3552 out.push_back(endIndices[i]);
3553 }
3554
3555 return out;
3556}
3557
3558vector<float> Scantable::getTsysSpectrum( int whichrow ) const
3559{
3560 Vector<Float> tsys( tsysCol_(whichrow) ) ;
3561 vector<float> stlTsys ;
3562 tsys.tovector( stlTsys ) ;
3563 return stlTsys ;
3564}
3565
3566
3567}
3568//namespace asap
Note: See TracBrowser for help on using the repository browser.