source: trunk/src/Scantable.cpp@ 2342

Last change on this file since 2342 was 2332, checked in by Takeshi Nakazato, 13 years ago

New Development: No

JIRA Issue: No

Ready for Test: No

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

###THIS SHOULD NOT BE INCLUDED RELEASE 3.3###
Take care of the data that the version number is 3 while the definition
of MOLECULE table is along version 4. For such data, just increment
version number.


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