Line data Source code
1 : #include <algorithm>
2 : #include <cstdlib>
3 : #include <iostream>
4 : #include <map>
5 : #include <string>
6 : #include <sys/time.h>
7 : #include <vector>
8 :
9 : #include <casacore/casa/Arrays/ArrayMath.h>
10 : #include <casacore/casa/Arrays/Vector.h>
11 : #include <casacore/casa/Containers/Allocator.h>
12 : #include <casacore/casa/Containers/Block.h>
13 : #include <casacore/casa/Logging/LogIO.h>
14 : #include <casacore/casa/Logging/LogOrigin.h>
15 : #include <casacore/casa/Quanta/MVTime.h>
16 : #include <casacore/casa/Utilities/Assert.h>
17 : #include <casacore/casa/Utilities/GenSort.h>
18 : #include <casacore/ms/MeasurementSets/MSSpectralWindow.h>
19 : #include <casacore/ms/MSSel/MSSelection.h>
20 : #include <casacore/ms/MSSel/MSSelectionTools.h>
21 : #include <msvis/MSVis/VisibilityIterator2.h>
22 : #include <msvis/MSVis/VisSetUtil.h>
23 : #include <casacore/scimath/Fitting/GenericL2Fit.h>
24 : #include <casacore/scimath/Fitting/NonLinearFitLM.h>
25 : #include <casacore/scimath/Functionals/CompiledFunction.h>
26 : #include <casacore/scimath/Functionals/CompoundFunction.h>
27 : #include <casacore/scimath/Functionals/Function.h>
28 : #include <casacore/scimath/Functionals/Gaussian1D.h>
29 : #include <casacore/scimath/Functionals/Lorentzian1D.h>
30 : #include <casacore/scimath/Mathematics/Convolver.h>
31 : #include <casacore/scimath/Mathematics/VectorKernel.h>
32 : #include <singledish/SingleDish/BaselineTable.h>
33 : #include <singledish/SingleDish/BLParameterParser.h>
34 : #include <singledish/SingleDish/LineFinder.h>
35 : #include <singledish/SingleDish/SingleDishMS.h>
36 : #include <stdcasa/StdCasa/CasacSupport.h>
37 : #include <casacore/tables/Tables/ScalarColumn.h>
38 : #include <casa_sakura/SakuraAlignedArray.h>
39 :
40 : // for importasap and importnro
41 : #include <singledishfiller/Filler/NRO2MSReader.h>
42 : #include <singledishfiller/Filler/Scantable2MSReader.h>
43 : #include <singledishfiller/Filler/SingleDishMSFiller.h>
44 :
45 : #define _ORIGIN LogOrigin("SingleDishMS", __func__, WHERE)
46 :
47 : namespace {
48 : // Max number of rows to get in each iteration
49 : constexpr casacore::Int kNRowBlocking = 1000;
50 : // Sinusoid
51 : constexpr int SinusoidWaveNumber_kUpperLimit = -999;
52 : // Weight
53 : constexpr size_t WeightIndex_kStddev = 0;
54 : constexpr size_t WeightIndex_kRms = 1;
55 : constexpr size_t WeightIndex_kNum = 2;
56 :
57 367 : double gettimeofday_sec() {
58 : struct timeval tv;
59 367 : gettimeofday(&tv, NULL);
60 367 : return tv.tv_sec + (double) tv.tv_usec * 1.0e-6;
61 : }
62 :
63 : using casa::vi::VisBuffer2;
64 : using casacore::Matrix;
65 : using casacore::Cube;
66 : using casacore::Float;
67 : using casacore::Complex;
68 : using casacore::AipsError;
69 :
70 : template<class CUBE_ACCESSOR>
71 : struct DataAccessorInterface {
72 962 : static void GetCube(VisBuffer2 const &vb, Cube<Float> &cube) {
73 962 : real(cube, CUBE_ACCESSOR::GetVisCube(vb));
74 962 : }
75 : static void GetSlice(VisBuffer2 const &vb, size_t const iPol,
76 : Matrix<Float> &cubeSlice) {
77 : real(cubeSlice, CUBE_ACCESSOR::GetVisCube(vb).yzPlane(iPol));
78 : }
79 : };
80 :
81 : struct DataAccessor: public DataAccessorInterface<DataAccessor> {
82 942 : static Cube<Complex> GetVisCube(VisBuffer2 const &vb) {
83 942 : return vb.visCube();
84 : }
85 : };
86 :
87 : struct CorrectedDataAccessor: public DataAccessorInterface<CorrectedDataAccessor> {
88 20 : static Cube<Complex> GetVisCube(VisBuffer2 const &vb) {
89 20 : return vb.visCubeCorrected();
90 : }
91 : };
92 :
93 : struct FloatDataAccessor {
94 987 : static void GetCube(VisBuffer2 const &vb, Cube<Float> &cube) {
95 987 : cube = GetVisCube(vb);
96 987 : }
97 : static void GetSlice(VisBuffer2 const &vb, size_t const iPol,
98 : Matrix<Float> &cubeSlice) {
99 : cubeSlice = GetVisCube(vb).yzPlane(iPol);
100 : }
101 : private:
102 987 : static Cube<Float> GetVisCube(VisBuffer2 const &vb) {
103 987 : return vb.visCubeFloat();
104 : }
105 : };
106 :
107 942 : inline void GetCubeFromData(VisBuffer2 const &vb, Cube<Float> &cube) {
108 942 : DataAccessor::GetCube(vb, cube);
109 942 : }
110 :
111 20 : inline void GetCubeFromCorrected(VisBuffer2 const &vb, Cube<Float> &cube) {
112 20 : CorrectedDataAccessor::GetCube(vb, cube);
113 20 : }
114 :
115 987 : inline void GetCubeFromFloat(VisBuffer2 const &vb, Cube<Float> &cube) {
116 987 : FloatDataAccessor::GetCube(vb, cube);
117 987 : }
118 :
119 0 : inline void GetCubeDefault(VisBuffer2 const& /*vb*/, Cube<Float>& /*cube*/) {
120 0 : throw AipsError("Data accessor for VB2 is not properly configured.");
121 : }
122 :
123 4708 : inline void compute_weight(size_t const num_data,
124 : float const data[/*num_data*/],
125 : bool const mask[/*num_data*/],
126 : std::vector<float>& weight) {
127 14124 : for (size_t i = 0; i < WeightIndex_kNum; ++i) {
128 9416 : weight[i] = 0.0;
129 : }
130 :
131 4708 : int num_data_effective = 0;
132 4708 : double sum = 0.0;
133 4708 : double sum_sq = 0.0;
134 19492964 : for (size_t i = 0; i < num_data; ++i) {
135 19488256 : if (mask[i]) {
136 19406262 : num_data_effective++;
137 19406262 : sum += data[i];
138 19406262 : sum_sq += data[i] * data[i];
139 : }
140 : }
141 :
142 4708 : if (num_data_effective > 0) {
143 4708 : double factor = 1.0 / static_cast<double>(num_data_effective);
144 4708 : double mean = sum * factor;
145 4708 : double mean_sq = sum_sq * factor;
146 :
147 9416 : std::vector<double> variance(WeightIndex_kNum);
148 4708 : variance[WeightIndex_kStddev] = mean_sq - mean * mean;
149 4708 : variance[WeightIndex_kRms] = mean_sq;
150 :
151 9416 : auto do_compute_weight = [&](size_t idx) {
152 9416 : if (variance[idx] > 0.0) {
153 9416 : weight[idx] = static_cast<float>(1.0 / variance[idx]);
154 : } else {
155 0 : LogIO os(_ORIGIN);
156 0 : os << "Weight set to 0 for a bad data." << LogIO::WARN;
157 : }
158 9416 : };
159 :
160 4708 : do_compute_weight(WeightIndex_kStddev);
161 4708 : do_compute_weight(WeightIndex_kRms);
162 : }
163 4708 : }
164 :
165 : } // anonymous namespace
166 :
167 : using namespace casacore;
168 : namespace casa {
169 :
170 0 : SingleDishMS::SingleDishMS() : msname_(""), sdh_(0) {
171 0 : initialize();
172 0 : }
173 :
174 754 : SingleDishMS::SingleDishMS(string const& ms_name) : msname_(ms_name), sdh_(0) {
175 2262 : LogIO os(_ORIGIN);
176 754 : initialize();
177 754 : }
178 :
179 754 : SingleDishMS::~SingleDishMS() {
180 754 : if (sdh_) {
181 12 : delete sdh_;
182 12 : sdh_ = 0;
183 : }
184 754 : msname_ = "";
185 754 : }
186 :
187 1490 : void SingleDishMS::initialize() {
188 1490 : in_column_ = MS::UNDEFINED_COLUMN;
189 : // out_column_ = MS::UNDEFINED_COLUMN;
190 1490 : doSmoothing_ = false;
191 1490 : doAtmCor_ = false;
192 1490 : visCubeAccessor_ = GetCubeDefault;
193 1490 : }
194 :
195 0 : bool SingleDishMS::close() {
196 0 : LogIO os(_ORIGIN);
197 0 : os << "Detaching from SingleDishMS" << LogIO::POST;
198 :
199 0 : if (sdh_) {
200 0 : delete sdh_;
201 0 : sdh_ = 0;
202 : }
203 0 : msname_ = "";
204 :
205 0 : return true;
206 : }
207 :
208 : ////////////////////////////////////////////////////////////////////////
209 : ///// Common utility functions
210 : ////////////////////////////////////////////////////////////////////////
211 753 : void SingleDishMS::setSelection(Record const &selection, bool const verbose) {
212 2259 : LogIO os(_ORIGIN);
213 753 : if (!selection_.empty()) // selection is set before
214 0 : os << "Discard old selection and setting new one." << LogIO::POST;
215 753 : if (selection.empty()) // empty selection is passed
216 0 : os << "Resetting selection." << LogIO::POST;
217 :
218 753 : selection_ = selection;
219 : // Verbose output
220 753 : bool any_selection(false);
221 753 : if (verbose && !selection_.empty()) {
222 1506 : String timeExpr(""), antennaExpr(""), fieldExpr(""), spwExpr(""),
223 1506 : uvDistExpr(""), taQLExpr(""), polnExpr(""), scanExpr(""), arrayExpr(""),
224 1506 : obsExpr(""), intentExpr("");
225 753 : timeExpr = get_field_as_casa_string(selection_, "timerange");
226 753 : antennaExpr = get_field_as_casa_string(selection_, "antenna");
227 753 : fieldExpr = get_field_as_casa_string(selection_, "field");
228 753 : spwExpr = get_field_as_casa_string(selection_, "spw");
229 753 : uvDistExpr = get_field_as_casa_string(selection_, "uvdist");
230 753 : taQLExpr = get_field_as_casa_string(selection_, "taql");
231 753 : polnExpr = get_field_as_casa_string(selection_, "correlation");
232 753 : scanExpr = get_field_as_casa_string(selection_, "scan");
233 753 : arrayExpr = get_field_as_casa_string(selection_, "array");
234 753 : intentExpr = get_field_as_casa_string(selection_, "intent");
235 753 : obsExpr = get_field_as_casa_string(selection_, "observation");
236 : // selection Summary
237 753 : os << "[Selection Summary]" << LogIO::POST;
238 753 : if (obsExpr != "") {
239 0 : any_selection = true;
240 0 : os << "- Observation: " << obsExpr << LogIO::POST;
241 : }
242 753 : if (antennaExpr != "") {
243 7 : any_selection = true;
244 7 : os << "- Antenna: " << antennaExpr << LogIO::POST;
245 : }
246 753 : if (fieldExpr != "") {
247 8 : any_selection = true;
248 8 : os << "- Field: " << fieldExpr << LogIO::POST;
249 : }
250 753 : if (spwExpr != "") {
251 578 : any_selection = true;
252 578 : os << "- SPW: " << spwExpr << LogIO::POST;
253 : }
254 753 : if (polnExpr != "") {
255 43 : any_selection = true;
256 43 : os << "- Pol: " << polnExpr << LogIO::POST;
257 : }
258 753 : if (scanExpr != "") {
259 9 : any_selection = true;
260 9 : os << "- Scan: " << scanExpr << LogIO::POST;
261 : }
262 753 : if (timeExpr != "") {
263 6 : any_selection = true;
264 6 : os << "- Time: " << timeExpr << LogIO::POST;
265 : }
266 753 : if (intentExpr != "") {
267 104 : any_selection = true;
268 104 : os << "- Intent: " << intentExpr << LogIO::POST;
269 : }
270 753 : if (arrayExpr != "") {
271 0 : any_selection = true;
272 0 : os << "- Array: " << arrayExpr << LogIO::POST;
273 : }
274 753 : if (uvDistExpr != "") {
275 0 : any_selection = true;
276 0 : os << "- UVDist: " << uvDistExpr << LogIO::POST;
277 : }
278 753 : if (taQLExpr != "") {
279 1 : any_selection = true;
280 1 : os << "- TaQL: " << taQLExpr << LogIO::POST;
281 : }
282 : {// reindex
283 : Int ifield;
284 753 : ifield = selection_.fieldNumber(String("reindex"));
285 753 : if (ifield > -1) {
286 753 : Bool reindex = selection_.asBool(ifield);
287 753 : os << "- Reindex: " << (reindex ? "ON" : "OFF" ) << LogIO::POST;
288 : }
289 : }
290 753 : if (!any_selection)
291 155 : os << "No valid selection parameter is set." << LogIO::WARN;
292 : }
293 753 : }
294 :
295 10 : void SingleDishMS::setAverage(Record const &average, bool const verbose) {
296 20 : LogIO os(_ORIGIN);
297 10 : if (!average_.empty()) // average is set before
298 0 : os << "Discard old average and setting new one." << LogIO::POST;
299 10 : if (average.empty()) // empty average is passed
300 0 : os << "Resetting average." << LogIO::POST;
301 :
302 10 : average_ = average;
303 :
304 10 : if (verbose && !average_.empty()) {
305 20 : LogIO os(_ORIGIN);
306 : Int ifield;
307 10 : ifield = average_.fieldNumber(String("timeaverage"));
308 10 : os << "[Averaging Settings]" << LogIO::POST;
309 10 : if (ifield < 0 || !average_.asBool(ifield)) {
310 0 : os << "No averaging will be done" << LogIO::POST;
311 0 : return;
312 : }
313 :
314 20 : String timebinExpr(""), timespanExpr(""), tweightExpr("");
315 10 : timebinExpr = get_field_as_casa_string(average_, "timebin");
316 10 : timespanExpr = get_field_as_casa_string(average_, "timespan");
317 10 : tweightExpr = get_field_as_casa_string(average_, "tweight");
318 10 : if (timebinExpr != "") {
319 10 : os << "- Time bin: " << timebinExpr << LogIO::POST;
320 : }
321 10 : if (timespanExpr != "") {
322 8 : os << "- Time span: " << timespanExpr << LogIO::POST;
323 : }
324 10 : if (tweightExpr != "") {
325 0 : os << "- Averaging weight: " << tweightExpr << LogIO::POST;
326 : }
327 :
328 : }
329 : }
330 :
331 4 : void SingleDishMS::setPolAverage(Record const &average, bool const verbose) {
332 8 : LogIO os(_ORIGIN);
333 4 : if (!pol_average_.empty()) // polaverage is set before
334 0 : os << "Discard old average and setting new one." << LogIO::POST;
335 4 : if (average.empty()) // empty polaverage is passed
336 0 : os << "Resetting average." << LogIO::POST;
337 :
338 4 : pol_average_ = average;
339 :
340 4 : if (verbose && !pol_average_.empty()) {
341 8 : LogIO os(_ORIGIN);
342 : Int ifield;
343 4 : ifield = pol_average_.fieldNumber(String("polaverage"));
344 4 : os << "[Polarization Averaging Settings]" << LogIO::POST;
345 4 : if (ifield < 0 || !pol_average_.asBool(ifield)) {
346 0 : os << "No polarization averaging will be done" << LogIO::POST;
347 0 : return;
348 : }
349 :
350 8 : String polAverageModeExpr("");
351 4 : polAverageModeExpr = get_field_as_casa_string(pol_average_, "polaveragemode");
352 4 : if (polAverageModeExpr != "") {
353 4 : os << "- Mode: " << polAverageModeExpr << LogIO::POST;
354 : }
355 : }
356 : }
357 :
358 9813 : String SingleDishMS::get_field_as_casa_string(Record const &in_data,
359 : string const &field_name) {
360 : Int ifield;
361 9813 : ifield = in_data.fieldNumber(String(field_name));
362 9813 : if (ifield > -1)
363 1358 : return in_data.asString(ifield);
364 8455 : return "";
365 : }
366 :
367 160 : bool SingleDishMS::prepare_for_process(string const &in_column_name,
368 : string const &out_ms_name) {
369 : // Sort by single dish default
370 160 : return prepare_for_process(in_column_name, out_ms_name, Block<Int>(), true);
371 : }
372 :
373 748 : bool SingleDishMS::prepare_for_process(string const &in_column_name,
374 : string const &out_ms_name,
375 : Block<Int> const &sortColumns,
376 : bool const addDefaultSortCols) {
377 2244 : LogIO os(_ORIGIN);
378 748 : AlwaysAssert(msname_ != "", AipsError);
379 : // define a column to read data from
380 1496 : string in_column_name_lower = in_column_name;
381 : std::transform(
382 : in_column_name_lower.begin(),
383 : in_column_name_lower.end(),
384 : in_column_name_lower.begin(),
385 5528 : [](unsigned char c) {return std::tolower(c);}
386 748 : );
387 748 : if (in_column_name_lower == "float_data") {
388 396 : in_column_ = MS::FLOAT_DATA;
389 396 : visCubeAccessor_ = GetCubeFromFloat;
390 352 : } else if (in_column_name_lower == "corrected") {
391 32 : in_column_ = MS::CORRECTED_DATA;
392 32 : visCubeAccessor_ = GetCubeFromCorrected;
393 320 : } else if (in_column_name_lower == "data") {
394 320 : in_column_ = MS::DATA;
395 320 : visCubeAccessor_ = GetCubeFromData;
396 : } else {
397 0 : throw(AipsError("Invalid data column name"));
398 : }
399 : // destroy SDMSManager
400 748 : if (sdh_)
401 0 : delete sdh_;
402 : // Configure record
403 1496 : Record configure_param(selection_);
404 748 : format_selection(configure_param);
405 748 : configure_param.define("inputms", msname_);
406 748 : configure_param.define("outputms", out_ms_name);
407 1496 : String in_name(in_column_name);
408 748 : in_name.upcase();
409 748 : configure_param.define("datacolumn", in_name);
410 : // merge averaging parameters
411 748 : configure_param.merge(average_);
412 : // The other available keys
413 : // - buffermode, realmodelcol, usewtspectrum, tileshape,
414 : // - chanaverage, chanbin, useweights,
415 : // - ddistart, hanning
416 : // - regridms, phasecenter, restfreq, outframe, interpolation, nspw,
417 : // - mode, nchan, start, width, veltype,
418 : // - timeaverage, timebin, timespan, maxuvwdistance
419 :
420 : // smoothing
421 748 : configure_param.define("smoothFourier", doSmoothing_);
422 :
423 : // merge polarization averaging parameters
424 748 : configure_param.merge(pol_average_);
425 :
426 : // offline ATM correction
427 748 : configure_param.define("atmCor", doAtmCor_);
428 748 : configure_param.merge(atmCorConfig_);
429 :
430 : // Generate SDMSManager
431 748 : sdh_ = new SDMSManager();
432 :
433 : // Configure SDMSManager
434 748 : sdh_->configure(configure_param);
435 :
436 1496 : ostringstream oss;
437 748 : configure_param.print(oss);
438 753 : String str(oss.str());
439 748 : os << LogIO::DEBUG1 << " Configuration Record " << LogIO::POST;
440 748 : os << LogIO::DEBUG1 << str << LogIO::POST;
441 : // Open the MS and select data
442 748 : sdh_->open();
443 743 : sdh_->getOutputMs()->flush();
444 : // set large timebin if not averaging
445 : Double timeBin;
446 743 : int exists = configure_param.fieldNumber("timebin");
447 743 : if (exists < 0) {
448 : // Avoid TIME col being added to sort columns in MSIter::construct.
449 : // TIME is automatically added to sort column when
450 : // timebin is not 0, even if addDefaultSortCols=false.
451 733 : timeBin = 0.0;
452 : } else {
453 20 : String timebin_string;
454 10 : configure_param.get(exists, timebin_string);
455 10 : timeBin = casaQuantity(timebin_string).get("s").getValue();
456 :
457 : Int ifield;
458 10 : ifield = configure_param.fieldNumber(String("timeaverage"));
459 10 : Bool average_time = ifield < 0 ? false : configure_param.asBool(ifield);
460 10 : if (timeBin < 0 || (average_time && timeBin == 0.0))
461 0 : throw(AipsError("time bin should be positive"));
462 : }
463 : // set sort column
464 743 : sdh_->setSortColumns(sortColumns, addDefaultSortCols, timeBin);
465 : // Set up the Data Handler
466 743 : sdh_->setup();
467 1486 : return true;
468 : }
469 :
470 736 : void SingleDishMS::finalize_process() {
471 736 : initialize();
472 736 : if (sdh_) {
473 736 : sdh_->close();
474 736 : delete sdh_;
475 736 : sdh_ = 0;
476 : }
477 736 : }
478 :
479 748 : void SingleDishMS::format_selection(Record &selection) {
480 : // At this moment sdh_ is not supposed to be generated yet.
481 2244 : LogIO os(_ORIGIN);
482 : // format spw
483 2244 : String const spwSel(get_field_as_casa_string(selection, "spw"));
484 923 : selection.define("spw", spwSel == "" ? "*" : spwSel);
485 :
486 : // Select only auto-correlation
487 1496 : String autoCorrSel("");
488 : os << "Formatting antenna selection to select only auto-correlation"
489 748 : << LogIO::POST;
490 1496 : String const antennaSel(get_field_as_casa_string(selection, "antenna"));
491 : os << LogIO::DEBUG1 << "Input antenna expression = " << antennaSel
492 748 : << LogIO::POST;
493 748 : if (antennaSel == "") { //Antenna selection is NOT set
494 741 : autoCorrSel = String("*&&&");
495 : } else { //User defined antenna selection
496 14 : MeasurementSet MSobj = MeasurementSet(msname_);
497 7 : MeasurementSet* theMS = &MSobj;
498 14 : MSSelection theSelection;
499 7 : theSelection.setAntennaExpr(antennaSel);
500 14 : TableExprNode exprNode = theSelection.toTableExprNode(theMS);
501 14 : Vector<Int> ant1Vec = theSelection.getAntenna1List();
502 : os << LogIO::DEBUG1 << ant1Vec.nelements()
503 7 : << " antenna(s) are selected. ID = ";
504 14 : for (uInt i = 0; i < ant1Vec.nelements(); ++i) {
505 7 : os << ant1Vec[i] << ", ";
506 7 : if (autoCorrSel != "")
507 0 : autoCorrSel += ";";
508 7 : autoCorrSel += String::toString(ant1Vec[i]) + "&&&";
509 : }
510 7 : os << LogIO::POST;
511 : }
512 : os << LogIO::DEBUG1 << "Auto-correlation selection string: " << autoCorrSel
513 748 : << LogIO::POST;
514 :
515 748 : selection.define("antenna", autoCorrSel);
516 748 : }
517 :
518 1949 : void SingleDishMS::get_data_cube_float(vi::VisBuffer2 const &vb,
519 : Cube<Float> &data_cube) {
520 : // if (in_column_ == MS::FLOAT_DATA) {
521 : // data_cube = vb.visCubeFloat();
522 : // } else { //need to convert Complex cube to Float
523 : // Cube<Complex> cdata_cube(data_cube.shape());
524 : // if (in_column_ == MS::DATA) {
525 : // cdata_cube = vb.visCube();
526 : // } else { //MS::CORRECTED_DATA
527 : // cdata_cube = vb.visCubeCorrected();
528 : // }
529 : // // convert Complext to Float
530 : // convertArrayC2F(data_cube, cdata_cube);
531 : // }
532 1949 : (*visCubeAccessor_)(vb, data_cube);
533 1949 : }
534 :
535 0 : void SingleDishMS::convertArrayC2F(Array<Float> &to,
536 : Array<Complex> const &from) {
537 0 : if (to.nelements() == 0 && from.nelements() == 0) {
538 0 : return;
539 : }
540 0 : if (to.shape() != from.shape()) {
541 0 : throw(ArrayConformanceError("Array shape differs"));
542 : }
543 0 : Array<Complex>::const_iterator endFrom = from.end();
544 0 : Array<Complex>::const_iterator iterFrom = from.begin();
545 0 : for (Array<Float>::iterator iterTo = to.begin(); iterFrom != endFrom; ++iterFrom, ++iterTo) {
546 0 : *iterTo = iterFrom->real();
547 : }
548 : }
549 :
550 39 : std::vector<string> SingleDishMS::split_string(string const &s, char delim) {
551 39 : std::vector<string> elems;
552 78 : string item;
553 86 : for (size_t i = 0; i < s.size(); ++i) {
554 47 : char ch = s.at(i);
555 47 : if (ch == delim) {
556 4 : if (!item.empty()) {
557 4 : elems.push_back(item);
558 : }
559 4 : item.clear();
560 : } else {
561 43 : item += ch;
562 : }
563 : }
564 39 : if (!item.empty()) {
565 39 : elems.push_back(item);
566 : }
567 78 : return elems;
568 : }
569 :
570 100 : bool SingleDishMS::file_exists(string const &filename) {
571 : FILE *fp;
572 100 : if ((fp = fopen(filename.c_str(), "r")) == NULL)
573 100 : return false;
574 0 : fclose(fp);
575 0 : return true;
576 : }
577 :
578 559 : void SingleDishMS::parse_spw(string const &in_spw,
579 : Vector<Int> &rec_spw,
580 : Matrix<Int> &rec_chan,
581 : Vector<size_t> &nchan,
582 : Vector<Vector<Bool> > &mask,
583 : Vector<bool> &nchan_set) {
584 1118 : Record selrec = sdh_->getSelRec(in_spw);
585 559 : rec_spw = selrec.asArrayInt("spw");
586 559 : rec_chan = selrec.asArrayInt("channel");
587 559 : nchan.resize(rec_spw.nelements());
588 559 : mask.resize(rec_spw.nelements());
589 559 : nchan_set.resize(rec_spw.nelements());
590 2569 : for (size_t i = 0; i < nchan_set.nelements(); ++i) {
591 2010 : nchan_set(i) = false;
592 : }
593 559 : }
594 :
595 1955 : void SingleDishMS::get_nchan_and_mask(Vector<Int> const &rec_spw,
596 : Vector<Int> const &data_spw,
597 : Matrix<Int> const &rec_chan,
598 : size_t const num_chan,
599 : Vector<size_t> &nchan,
600 : Vector<Vector<Bool> > &mask,
601 : Vector<bool> &nchan_set,
602 : bool &new_nchan) {
603 1955 : new_nchan = false;
604 7300 : for (size_t i = 0; i < rec_spw.nelements(); ++i) {
605 : //get nchan by spwid and set to nchan[]
606 748521 : for (size_t j = 0; j < data_spw.nelements(); ++j) {
607 744342 : if ((!nchan_set(i)) && (data_spw(j) == rec_spw(i))) {
608 1166 : bool found = false;
609 5418 : for (size_t k = 0; k < nchan.nelements(); ++k) {
610 4252 : if (!nchan_set(k))
611 3146 : continue;
612 1106 : if (nchan(k) == num_chan)
613 1106 : found = true;
614 : }
615 1166 : if (!found) {
616 559 : new_nchan = true;
617 : }
618 1166 : nchan(i) = num_chan;
619 1166 : nchan_set(i) = true;
620 1166 : break;
621 : }
622 : }
623 5345 : if (!nchan_set(i))
624 2284 : continue;
625 3061 : mask(i).resize(nchan(i));
626 : // generate mask
627 3061 : get_mask_from_rec(rec_spw(i), rec_chan, mask(i), true);
628 : }
629 1955 : }
630 :
631 3138 : void SingleDishMS::get_mask_from_rec(Int spwid,
632 : Matrix<Int> const &rec_chan,
633 : Vector<Bool> &mask,
634 : bool initialize) {
635 3138 : if (initialize) {
636 17243209 : for (size_t j = 0; j < mask.nelements(); ++j) {
637 17240148 : mask(j) = false;
638 : }
639 : }
640 : //construct a list of (start, end, stride, start, end, stride, ...)
641 : //from rec_chan for the spwid
642 6276 : std::vector<uInt> edge;
643 3138 : edge.clear();
644 13696 : for (size_t j = 0; j < rec_chan.nrow(); ++j) {
645 10558 : if (rec_chan.row(j)(0) == spwid) {
646 3696 : edge.push_back(rec_chan.row(j)(1));
647 3696 : edge.push_back(rec_chan.row(j)(2));
648 3696 : edge.push_back(rec_chan.row(j)(3));
649 : }
650 : }
651 : //generate mask
652 6834 : for (size_t j = 0; j < edge.size()-2; j += 3) {
653 17464145 : for (size_t k = edge[j]; k <= edge[j + 1] && k < mask.size(); k += edge[j + 2]) {
654 17460449 : mask(k) = true;
655 : }
656 : }
657 3138 : }
658 :
659 798635 : void SingleDishMS::get_masklist_from_mask(size_t const num_chan,
660 : bool const *mask,
661 : Vector<uInt> &masklist) {
662 798635 : size_t const max_num_masklist = num_chan + 1;
663 798635 : masklist.resize(max_num_masklist); // clear
664 798635 : uInt last_idx = num_chan - 1;
665 798635 : uInt num_masklist = 0;
666 2404532 : auto append = [&](uInt i){
667 2404532 : masklist[num_masklist] = i;
668 2404532 : num_masklist++;
669 3203167 : };
670 :
671 798635 : if (mask[0]) {
672 267551 : append(0);
673 : }
674 66435221 : for (uInt i = 1; i < last_idx; ++i) {
675 65636586 : if (!mask[i]) continue;
676 :
677 : // The following if-statements must be judged independently.
678 : // Don't put them together as a single statement by connecting with '||'.
679 55538636 : if (!mask[i - 1]) {
680 934697 : append(i);
681 : }
682 55538636 : if (!mask[i + 1]) {
683 934718 : append(i);
684 : }
685 : }
686 798635 : if (mask[last_idx]) {
687 267548 : if ((1 <= last_idx) && (!mask[last_idx - 1])) {
688 18 : append(last_idx);
689 : }
690 267548 : append(last_idx);
691 : }
692 798635 : masklist.resize(num_masklist, true);
693 798635 : }
694 :
695 351 : void SingleDishMS::get_baseline_context(size_t const bltype,
696 : uint16_t order,
697 : size_t num_chan,
698 : Vector<size_t> const &nchan,
699 : Vector<bool> const &nchan_set,
700 : Vector<size_t> &ctx_indices,
701 : std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> &bl_contexts) {
702 351 : size_t idx = 0;
703 351 : bool found = false;
704 445 : for (size_t i = 0; i < nchan.nelements(); ++i) {
705 445 : if ((nchan_set[i])&&(nchan[i] == num_chan)) {
706 351 : idx = bl_contexts.size();
707 351 : found = true;
708 351 : break;
709 : }
710 : }
711 351 : if (found) {
712 1426 : for (size_t i = 0; i < nchan.nelements(); ++i) {
713 1075 : if ((nchan_set[i])&&(nchan[i] == num_chan)) {
714 351 : ctx_indices[i] = idx;
715 : }
716 : }
717 :
718 : LIBSAKURA_SYMBOL(LSQFitContextFloat) *context;
719 351 : LIBSAKURA_SYMBOL(Status) status = LIBSAKURA_SYMBOL(Status_kNG);
720 351 : if ((bltype == BaselineType_kPolynomial)||(bltype == BaselineType_kChebyshev)) {
721 241 : status = LIBSAKURA_SYMBOL(CreateLSQFitContextPolynomialFloat)(static_cast<LIBSAKURA_SYMBOL(LSQFitType)>(bltype),
722 : static_cast<uint16_t>(order),
723 : num_chan, &context);
724 110 : } else if (bltype == BaselineType_kCubicSpline) {
725 87 : status = LIBSAKURA_SYMBOL(CreateLSQFitContextCubicSplineFloat)(static_cast<uint16_t>(order),
726 : num_chan, &context);
727 23 : } else if (bltype == BaselineType_kSinusoid) {
728 23 : status = LIBSAKURA_SYMBOL(CreateLSQFitContextSinusoidFloat)(static_cast<uint16_t>(order),
729 : num_chan, &context);
730 : }
731 351 : check_sakura_status("sakura_CreateLSQFitContextFloat", status);
732 351 : bl_contexts.push_back(context);
733 : }
734 351 : }
735 :
736 371 : void SingleDishMS::get_baseline_context(size_t const bltype,
737 : uint16_t order,
738 : size_t num_chan,
739 : size_t ispw,
740 : Vector<size_t> &ctx_indices,
741 : std::vector<size_t> & ctx_nchans,
742 : std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> &bl_contexts) {
743 371 : AlwaysAssert(bl_contexts.size() == ctx_nchans.size() || bl_contexts.size() == ctx_nchans.size()-1 , AipsError);
744 371 : size_t idx = 0;
745 371 : bool found = false;
746 371 : for (size_t i = 0; i < bl_contexts.size(); ++i) {
747 211 : if (ctx_nchans[i] == num_chan) {
748 211 : idx = i;
749 211 : found = true;
750 211 : break;
751 : }
752 : }
753 371 : if (found) {
754 : // contexts with the valid number of channels already exists.
755 : // just update idx to bl_contexts and return.
756 211 : ctx_indices[ispw] = idx;
757 211 : return;
758 : }
759 : // contexts with the number of channels is not yet in bl_contexts.
760 : // Need to create a new context.
761 160 : ctx_indices[ispw] = bl_contexts.size();
762 : LIBSAKURA_SYMBOL(LSQFitContextFloat) *context;
763 160 : LIBSAKURA_SYMBOL(Status) status = LIBSAKURA_SYMBOL(Status_kNG);
764 160 : if ((bltype == BaselineType_kPolynomial)||(bltype == BaselineType_kChebyshev)) {
765 110 : status = LIBSAKURA_SYMBOL(CreateLSQFitContextPolynomialFloat)(static_cast<LIBSAKURA_SYMBOL(LSQFitType)>(bltype),
766 : static_cast<uint16_t>(order),
767 : num_chan, &context);
768 50 : } else if (bltype == BaselineType_kCubicSpline) {
769 50 : status = LIBSAKURA_SYMBOL(CreateLSQFitContextCubicSplineFloat)(static_cast<uint16_t>(order),
770 : num_chan, &context);
771 0 : } else if (bltype == BaselineType_kSinusoid) {
772 0 : status = LIBSAKURA_SYMBOL(CreateLSQFitContextSinusoidFloat)(static_cast<uint16_t>(order),
773 : num_chan, &context);
774 : }
775 160 : check_sakura_status("sakura_CreateLSQFitContextFloat", status);
776 160 : bl_contexts.push_back(context);
777 160 : if (ctx_nchans.size() != bl_contexts.size()) {
778 57 : ctx_nchans.push_back(num_chan);
779 : }
780 160 : AlwaysAssert(bl_contexts.size() == ctx_nchans.size(), AipsError);
781 : }
782 :
783 640 : void SingleDishMS::destroy_baseline_contexts(std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> &bl_contexts) {
784 : LIBSAKURA_SYMBOL(Status) status;
785 1151 : for (size_t i = 0; i < bl_contexts.size(); ++i) {
786 511 : status = LIBSAKURA_SYMBOL(DestroyLSQFitContextFloat)(bl_contexts[i]);
787 511 : check_sakura_status("sakura_DestoyBaselineContextFloat", status);
788 : }
789 640 : }
790 :
791 2397975 : void SingleDishMS::check_sakura_status(string const &name,
792 : LIBSAKURA_SYMBOL(Status) const status) {
793 2397975 : if (status == LIBSAKURA_SYMBOL(Status_kOK)) return;
794 :
795 0 : ostringstream oss;
796 0 : oss << name << "() failed -- ";
797 0 : if (status == LIBSAKURA_SYMBOL(Status_kNG)) {
798 0 : oss << "NG";
799 0 : } else if (status == LIBSAKURA_SYMBOL(Status_kInvalidArgument)) {
800 0 : oss << "InvalidArgument";
801 0 : } else if (status == LIBSAKURA_SYMBOL(Status_kNoMemory)) {
802 0 : oss << "NoMemory";
803 0 : } else if (status == LIBSAKURA_SYMBOL(Status_kUnknownError)) {
804 0 : oss << "UnknownError";
805 : }
806 0 : throw(AipsError(oss.str()));
807 : }
808 :
809 198610 : void SingleDishMS::check_baseline_status(LIBSAKURA_SYMBOL(LSQFitStatus) const bl_status) {
810 198610 : if (bl_status != LIBSAKURA_SYMBOL(LSQFitStatus_kOK)) {
811 0 : throw(AipsError("baseline fitting isn't successful."));
812 : }
813 198610 : }
814 :
815 799515 : void SingleDishMS::get_spectrum_from_cube(Cube<Float> &data_cube,
816 : size_t const row,
817 : size_t const plane,
818 : size_t const num_data,
819 : float out_data[]) {
820 799515 : AlwaysAssert(static_cast<size_t>(data_cube.ncolumn()) == num_data, AipsError);
821 74645595 : for (size_t i = 0; i < num_data; ++i)
822 73846080 : out_data[i] = static_cast<float>(data_cube(plane, i, row));
823 799515 : }
824 :
825 799309 : void SingleDishMS::set_spectrum_to_cube(Cube<Float> &data_cube,
826 : size_t const row,
827 : size_t const plane,
828 : size_t const num_data,
829 : float *in_data) {
830 799309 : AlwaysAssert(static_cast<size_t>(data_cube.ncolumn()) == num_data, AipsError);
831 73554573 : for (size_t i = 0; i < num_data; ++i)
832 72755264 : data_cube(plane, i, row) = static_cast<Float>(in_data[i]);
833 799309 : }
834 :
835 37 : void SingleDishMS::get_weight_matrix(vi::VisBuffer2 const &vb,
836 : Matrix<Float> &weight_matrix) {
837 37 : weight_matrix = vb.weight();
838 37 : }
839 :
840 4708 : void SingleDishMS::set_weight_to_matrix(Matrix<Float> &weight_matrix,
841 : size_t const row,
842 : size_t const plane,
843 : float in_weight) {
844 4708 : weight_matrix(plane, row) = static_cast<Float>(in_weight);
845 4708 : }
846 :
847 1949 : void SingleDishMS::get_flag_cube(vi::VisBuffer2 const &vb,
848 : Cube<Bool> &flag_cube) {
849 1949 : flag_cube = vb.flagCube();
850 1949 : }
851 :
852 805597 : void SingleDishMS::get_flag_from_cube(Cube<Bool> &flag_cube,
853 : size_t const row,
854 : size_t const plane,
855 : size_t const num_flag,
856 : bool out_flag[]) {
857 805597 : AlwaysAssert(static_cast<size_t>(flag_cube.ncolumn()) == num_flag, AipsError);
858 100099485 : for (size_t i = 0; i < num_flag; ++i)
859 99293888 : out_flag[i] = static_cast<bool>(flag_cube(plane, i, row));
860 805597 : }
861 :
862 0 : void SingleDishMS::set_flag_to_cube(Cube<Bool> &flag_cube,
863 : size_t const row,
864 : size_t const plane,
865 : size_t const num_flag,
866 : bool *in_flag) {
867 0 : AlwaysAssert(static_cast<size_t>(flag_cube.ncolumn()) == num_flag, AipsError);
868 0 : for (size_t i = 0; i < num_flag; ++i)
869 0 : flag_cube(plane, i, row) = static_cast<Bool>(in_flag[i]);
870 0 : }
871 :
872 1676 : void SingleDishMS::flag_spectrum_in_cube(Cube<Bool> &flag_cube,
873 : size_t const row,
874 : size_t const plane) {
875 1676 : uInt const num_flag = flag_cube.ncolumn();
876 7415436 : for (uInt ichan = 0; ichan < num_flag; ++ichan)
877 7413760 : flag_cube(plane, ichan, row) = true;
878 1676 : }
879 :
880 805597 : bool SingleDishMS::allchannels_flagged(size_t const num_flag,
881 : bool const* flag) {
882 805597 : bool res = true;
883 20583521 : for (size_t i = 0; i < num_flag; ++i) {
884 20579112 : if (!flag[i]) {
885 801188 : res = false;
886 801188 : break;
887 : }
888 : }
889 805597 : return res;
890 : }
891 :
892 799558 : size_t SingleDishMS::NValidMask(size_t const num_mask, bool const* mask) {
893 799558 : std::size_t nvalid = 0;
894 : // the assertion lines had better be replaced with static_assert when c++11 is supported
895 799558 : AlwaysAssert(static_cast<std::size_t>(true) == 1, AipsError);
896 799558 : AlwaysAssert(static_cast<std::size_t>(false) == 0, AipsError);
897 75594630 : for (size_t i = 0; i < num_mask; ++i) {
898 74795072 : nvalid += static_cast<std::size_t>(mask[i]);
899 : }
900 799558 : return nvalid;
901 : }
902 :
903 499 : void SingleDishMS::split_bloutputname(string str) {
904 499 : char key = ',';
905 998 : vector<size_t> v;
906 14736 : for (size_t i = 0; i < str.size(); ++i) {
907 14237 : char target = str[i];
908 14237 : if (key == target) {
909 998 : v.push_back(i);
910 : }
911 : }
912 :
913 : //cout << "comma " << v.size() << endl;
914 : //cout << "v[1] " << v[1] << endl;
915 : //cout << "v.size()-1 " << v.size()-1 << endl;
916 : //cout << "v[1]+1 " << v[1]+1 << endl;
917 : //cout << "str.size()-v[1]-1 " << str.size()-v[1]-1 << endl;
918 : //cout << "str.substr(v[1]+1, str.size()-v[1]-1) " << str.substr(v[1]+1, str.size()-v[1]-1) << endl;
919 :
920 998 : string ss;
921 499 : bloutputname_csv.clear();
922 499 : bloutputname_text.clear();
923 499 : bloutputname_table.clear();
924 :
925 499 : if (0 != v[0]) {
926 222 : bloutputname_csv = str.substr(0, v[0]);
927 222 : ss = str.substr(0, v[0]);
928 : }
929 499 : if (v[0] + 1 != v[1]) {
930 100 : bloutputname_text = str.substr(v[0] + 1, v[1] - v[0] - 1);
931 : }
932 499 : if (v[1] != str.size() - 1) {
933 74 : bloutputname_table = str.substr(v[1] + 1, str.size() - v[1] - 1);
934 : }
935 499 : }
936 :
937 798635 : size_t SingleDishMS::get_num_coeff_bloutput(size_t const bltype,
938 : size_t order,
939 : size_t &num_coeff_max) {
940 798635 : size_t num_coeff = 0;
941 798635 : switch (bltype) {
942 401487 : case BaselineType_kPolynomial:
943 : case BaselineType_kChebyshev:
944 401487 : break;
945 198706 : case BaselineType_kCubicSpline:
946 198706 : num_coeff = order + 1;
947 198706 : break;
948 198442 : case BaselineType_kSinusoid:
949 198442 : break;
950 0 : default:
951 0 : throw(AipsError("Unsupported baseline type."));
952 : }
953 798635 : if (num_coeff_max < num_coeff) {
954 198642 : num_coeff_max = num_coeff;
955 : }
956 :
957 798635 : return num_coeff;
958 : }
959 :
960 228 : vector<int> SingleDishMS::string_to_list(string const &wn_str, char const delim) {
961 228 : vector<int> wn_list;
962 228 : wn_list.clear();
963 456 : vector<size_t> delim_position;
964 228 : delim_position.clear();
965 24941 : for (size_t i = 0; i < wn_str.size(); ++i) {
966 24713 : if (wn_str[i] == delim) {
967 5168 : delim_position.push_back(i);
968 : }
969 : }
970 228 : delim_position.push_back(wn_str.size());
971 228 : size_t start_position = 0;
972 5624 : for (size_t i = 0; i < delim_position.size(); ++i) {
973 5396 : size_t end_position = delim_position[i];
974 5396 : size_t length = end_position - start_position;
975 5396 : if (length > 0) {
976 5288 : wn_list.push_back(std::atoi(wn_str.substr(start_position, length).c_str()));
977 : }
978 5396 : start_position = end_position + 1;
979 : }
980 456 : return wn_list;
981 : }
982 :
983 110 : void SingleDishMS::get_effective_nwave(std::vector<int> const &addwn,
984 : std::vector<int> const &rejwn,
985 : int const wn_ulimit,
986 : std::vector<int> &effwn) {
987 110 : effwn.clear();
988 110 : if (addwn.size() < 1) {
989 0 : throw AipsError("addwn has no elements.");
990 : }
991 :
992 220 : auto up_to_nyquist_limit = [&](std::vector<int> const &v){ return ((v.size() == 2) && (v[1] == SinusoidWaveNumber_kUpperLimit)); };
993 266 : auto check_rejwn_add = [&](int const v){
994 266 : bool add = (0 <= v) && (v <= wn_ulimit); // check v in range
995 393 : for (size_t i = 0; i < rejwn.size(); ++i) {
996 148 : if (rejwn[i] == v) {
997 21 : add = false;
998 21 : break;
999 : }
1000 : }
1001 266 : if (add) {
1002 243 : effwn.push_back(v);
1003 : }
1004 376 : };
1005 :
1006 110 : if (up_to_nyquist_limit(addwn)) {
1007 8 : if (up_to_nyquist_limit(rejwn)) {
1008 0 : if (addwn[0] < rejwn[0]) {
1009 0 : for (int wn = addwn[0]; wn < rejwn[0]; ++wn) {
1010 0 : if ((0 <= wn) && (wn <= wn_ulimit)) {
1011 0 : effwn.push_back(wn);
1012 : }
1013 : }
1014 : } else {
1015 0 : throw(AipsError("No effective wave number given for sinusoidal fitting."));
1016 : }
1017 : } else {
1018 137 : for (int wn = addwn[0]; wn <= wn_ulimit; ++wn) {
1019 129 : check_rejwn_add(wn);
1020 : }
1021 : }
1022 : } else {
1023 102 : if (up_to_nyquist_limit(rejwn)) {
1024 1 : int maxwn = rejwn[0] - 1;
1025 1 : if (maxwn < 0) {
1026 0 : throw(AipsError("No effective wave number given for sinusoidal fitting."));
1027 : }
1028 3 : for (size_t i = 0; i < addwn.size(); ++i) {
1029 2 : if ((0 <= addwn[i]) && (addwn[i] <= maxwn)) {
1030 1 : effwn.push_back(addwn[i]);
1031 : }
1032 : }
1033 : } else {
1034 238 : for (size_t i = 0; i < addwn.size(); ++i) {
1035 137 : check_rejwn_add(addwn[i]);
1036 : }
1037 : }
1038 : }
1039 110 : if (effwn.size() == 0) {
1040 6 : throw(AipsError("No effective wave number given for sinusoidal fitting."));
1041 : }
1042 104 : }
1043 :
1044 198610 : void SingleDishMS::finalise_effective_nwave(std::vector<int> const &blparam_eff_base,
1045 : std::vector<int> const &blparam_exclude,
1046 : int const &blparam_upperlimit,
1047 : size_t const &num_chan,
1048 : float const *spec,
1049 : bool const *mask,
1050 : bool const &applyfft,
1051 : string const &fftmethod,
1052 : string const &fftthresh_str,
1053 : std::vector<size_t> &blparam_eff) {
1054 198610 : blparam_eff.resize(blparam_eff_base.size());
1055 198610 : copy(blparam_eff_base.begin(), blparam_eff_base.end(), blparam_eff.begin());
1056 :
1057 198610 : if (applyfft) {
1058 397128 : string fftthresh_attr;
1059 : float fftthresh_sigma;
1060 : int fftthresh_top;
1061 198564 : parse_fftthresh(fftthresh_str, fftthresh_attr, fftthresh_sigma, fftthresh_top);
1062 397128 : std::vector<int> blparam_fft;
1063 198564 : select_wavenumbers_via_fft(num_chan, spec, mask, fftmethod, fftthresh_attr,
1064 : fftthresh_sigma, fftthresh_top, blparam_upperlimit,
1065 : blparam_fft);
1066 198564 : merge_wavenumbers(blparam_eff_base, blparam_fft, blparam_exclude, blparam_eff);
1067 : }
1068 198610 : }
1069 :
1070 198564 : void SingleDishMS::parse_fftthresh(string const& fftthresh_str,
1071 : string& fftthresh_attr,
1072 : float& fftthresh_sigma,
1073 : int& fftthresh_top) {
1074 198564 : size_t idx_sigma = fftthresh_str.find("sigma");
1075 198564 : size_t idx_top = fftthresh_str.find("top");
1076 :
1077 198564 : if (idx_top == 0) {
1078 8 : std::istringstream is(fftthresh_str.substr(3));
1079 4 : is >> fftthresh_top;
1080 4 : fftthresh_attr = "top";
1081 198560 : } else if (idx_sigma == fftthresh_str.size() - 5) {
1082 12 : std::istringstream is(fftthresh_str.substr(0, fftthresh_str.size() - 5));
1083 6 : is >> fftthresh_sigma;
1084 6 : fftthresh_attr = "sigma";
1085 : } else {
1086 198554 : bool is_number = true;
1087 1588432 : for (size_t i = 0; i < fftthresh_str.size()-1; ++i) {
1088 1389878 : char ch = (fftthresh_str.substr(i, 1).c_str())[0];
1089 1389878 : if (!(isdigit(ch) || (fftthresh_str.substr(i, 1) == "."))) {
1090 0 : is_number = false;
1091 0 : break;
1092 : }
1093 : }
1094 198554 : if (is_number) {
1095 397108 : std::istringstream is(fftthresh_str);
1096 198554 : is >> fftthresh_sigma;
1097 198554 : fftthresh_attr = "sigma";
1098 : } else {
1099 0 : throw(AipsError("fftthresh has a wrong value"));
1100 : }
1101 : }
1102 198564 : }
1103 :
1104 198564 : void SingleDishMS::select_wavenumbers_via_fft(size_t const num_chan,
1105 : float const *spec,
1106 : bool const *mask,
1107 : string const &fftmethod,
1108 : string const &fftthresh_attr,
1109 : float const fftthresh_sigma,
1110 : int const fftthresh_top,
1111 : int const blparam_upperlimit,
1112 : std::vector<int> &blparam_fft) {
1113 198564 : blparam_fft.clear();
1114 397128 : std::vector<float> fourier_spec;
1115 198564 : if (fftmethod == "fft") {
1116 198564 : exec_fft(num_chan, spec, mask, false, true, fourier_spec);
1117 : } else {
1118 0 : throw AipsError("fftmethod must be 'fft' for now.");
1119 : }
1120 :
1121 198564 : int fourier_spec_size = static_cast<int>(fourier_spec.size());
1122 198564 : if (fftthresh_attr == "sigma") {
1123 198560 : float mean = 0.0;
1124 198560 : float mean2 = 0.0;
1125 6937408 : for (int i = 0; i < fourier_spec_size; ++i) {
1126 6738848 : mean += fourier_spec[i];
1127 6738848 : mean2 += fourier_spec[i] * fourier_spec[i];
1128 : }
1129 198560 : mean /= static_cast<float>(fourier_spec_size);
1130 198560 : mean2 /= static_cast<float>(fourier_spec_size);
1131 198560 : float thres = mean + fftthresh_sigma * float(sqrt(mean2 - mean * mean));
1132 :
1133 6937408 : for (int i = 0; i < fourier_spec_size; ++i) {
1134 6738848 : if ((i <= blparam_upperlimit)&&(thres <= fourier_spec[i])) {
1135 203720 : blparam_fft.push_back(i);
1136 : }
1137 : }
1138 4 : } else if (fftthresh_attr == "top") {
1139 4 : int i = 0;
1140 20 : while (i < fftthresh_top) {
1141 16 : float max = 0.0;
1142 16 : int max_idx = 0;
1143 65568 : for (int j = 0; j < fourier_spec_size; ++j) {
1144 65552 : if (max < fourier_spec[j]) {
1145 32 : max = fourier_spec[j];
1146 32 : max_idx = j;
1147 : }
1148 : }
1149 16 : fourier_spec[max_idx] = 0.0;
1150 16 : if (max_idx <= blparam_upperlimit) {
1151 16 : blparam_fft.push_back(max_idx);
1152 16 : ++i;
1153 : }
1154 : }
1155 : } else {
1156 0 : throw AipsError("fftthresh is wrong.");
1157 : }
1158 198564 : }
1159 :
1160 198564 : void SingleDishMS::exec_fft(size_t const num_chan,
1161 : float const *in_spec,
1162 : bool const *in_mask,
1163 : bool const get_real_imag,
1164 : bool const get_ampl_only,
1165 : std::vector<float> &fourier_spec) {
1166 397128 : Vector<Float> spec;
1167 198564 : interpolate_constant(static_cast<int>(num_chan), in_spec, in_mask, spec);
1168 :
1169 397128 : FFTServer<Float, Complex> ffts;
1170 397128 : Vector<Complex> fftres;
1171 198564 : ffts.fft0(fftres, spec);
1172 :
1173 198564 : float norm = static_cast<float>(2.0/static_cast<double>(num_chan));
1174 198564 : fourier_spec.clear();
1175 198564 : if (get_real_imag) {
1176 0 : for (size_t i = 0; i < fftres.size(); ++i) {
1177 0 : fourier_spec.push_back(real(fftres[i]) * norm);
1178 0 : fourier_spec.push_back(imag(fftres[i]) * norm);
1179 : }
1180 : } else {
1181 6953800 : for (size_t i = 0; i < fftres.size(); ++i) {
1182 6755236 : fourier_spec.push_back(abs(fftres[i]) * norm);
1183 6755236 : if (!get_ampl_only) fourier_spec.push_back(arg(fftres[i]));
1184 : }
1185 : }
1186 198564 : }
1187 :
1188 198564 : void SingleDishMS::interpolate_constant(int const num_chan,
1189 : float const *in_spec,
1190 : bool const *in_mask,
1191 : Vector<Float> &spec) {
1192 198564 : spec.resize(num_chan);
1193 13311908 : for (int i = 0; i < num_chan; ++i) {
1194 13113344 : spec[i] = in_spec[i];
1195 : }
1196 198564 : int idx_left = -1;
1197 198564 : int idx_right = -1;
1198 198564 : bool masked_region = false;
1199 :
1200 10730788 : for (int i = 0; i < num_chan; ++i) {
1201 10532224 : if (!in_mask[i]) {
1202 366080 : masked_region = true;
1203 366080 : idx_left = i;
1204 3079936 : while (i < num_chan) {
1205 2947200 : if (in_mask[i]) break;
1206 2713856 : idx_right = i;
1207 2713856 : ++i;
1208 : }
1209 : }
1210 :
1211 10532224 : if (masked_region) {
1212 : // execute interpolation as the following criteria:
1213 : // (1) for a masked region inside the spectrum, replace the spectral
1214 : // values with the mean of those at the two channels just outside
1215 : // the both edges of the masked region.
1216 : // (2) for a masked region at the spectral edge, replace the values
1217 : // with the one at the nearest non-masked channel.
1218 : // (ZOH, but bilateral)
1219 366080 : Float interp = 0.0;
1220 366080 : int idx_left_next = idx_left - 1;
1221 366080 : int idx_right_next = idx_right + 1;
1222 366080 : if (idx_left_next < 0) {
1223 132736 : if (idx_right_next < num_chan) {
1224 132736 : interp = in_spec[idx_right_next];
1225 : } else {
1226 0 : throw AipsError("Bad data. all channels are masked.");
1227 : }
1228 : } else {
1229 233344 : interp = in_spec[idx_left_next];
1230 233344 : if (idx_right_next < num_chan) {
1231 100608 : interp = (interp + in_spec[idx_right_next]) / 2.0;
1232 : }
1233 : }
1234 :
1235 366080 : if ((0 <= idx_left) && (idx_left < num_chan) && (0 <= idx_right) && (idx_right < num_chan)) {
1236 3079936 : for (int j = idx_left; j <= idx_right; ++j) {
1237 2713856 : spec[j] = interp;
1238 : }
1239 : }
1240 :
1241 366080 : masked_region = false;
1242 : }
1243 : }
1244 198564 : }
1245 :
1246 198564 : void SingleDishMS::merge_wavenumbers(std::vector<int> const &blparam_eff_base,
1247 : std::vector<int> const &blparam_fft,
1248 : std::vector<int> const &blparam_exclude,
1249 : std::vector<size_t> &blparam_eff) {
1250 402300 : for (size_t i = 0; i < blparam_fft.size(); ++i) {
1251 203736 : bool found = false;
1252 208944 : for (size_t j = 0; j < blparam_eff_base.size(); ++j) {
1253 203744 : if (blparam_eff_base[j] == blparam_fft[i]) {
1254 198536 : found = true;
1255 198536 : break;
1256 : }
1257 : }
1258 203736 : if (!found) { //new value to add
1259 : //but still need to check if it is to be excluded
1260 5200 : bool found_exclude = false;
1261 5200 : for (size_t j = 0; j < blparam_exclude.size(); ++j) {
1262 16 : if (blparam_exclude[j] == blparam_fft[i]) {
1263 16 : found_exclude = true;
1264 16 : break;
1265 : }
1266 : }
1267 5200 : if (!found_exclude) {
1268 5184 : blparam_eff.push_back(blparam_fft[i]);
1269 : }
1270 : }
1271 : }
1272 :
1273 198564 : if (1 < blparam_eff.size()) {
1274 928 : sort(blparam_eff.begin(), blparam_eff.end());
1275 928 : unique(blparam_eff.begin(), blparam_eff.end());
1276 : }
1277 198564 : }
1278 :
1279 : template<typename Func0, typename Func1, typename Func2, typename Func3>
1280 418 : void SingleDishMS::doSubtractBaseline(string const& in_column_name,
1281 : string const& out_ms_name,
1282 : string const& out_bloutput_name,
1283 : bool const& do_subtract,
1284 : string const& in_spw,
1285 : bool const& update_weight,
1286 : string const& sigma_value,
1287 : LIBSAKURA_SYMBOL(Status)& status,
1288 : std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> &bl_contexts,
1289 : size_t const bltype,
1290 : vector<int> const& blparam,
1291 : vector<int> const& blparam_exclude,
1292 : bool const& applyfft,
1293 : string const& fftmethod,
1294 : string const& fftthresh,
1295 : float const clip_threshold_sigma,
1296 : int const num_fitting_max,
1297 : bool const linefinding,
1298 : float const threshold,
1299 : int const avg_limit,
1300 : int const minwidth,
1301 : vector<int> const& edge,
1302 : Func0 func0,
1303 : Func1 func1,
1304 : Func2 func2,
1305 : Func3 func3,
1306 : LogIO os) {
1307 418 : os << LogIO::DEBUG2 << "Calling SingleDishMS::doSubtractBaseline " << LogIO::POST;
1308 :
1309 : // in_ms = out_ms
1310 : // in_column = [FLOAT_DATA|DATA|CORRECTED_DATA], out_column=new MS
1311 : // no iteration is necessary for the processing.
1312 : // procedure
1313 : // 1. iterate over MS
1314 : // 2. pick single spectrum from in_column (this is not actually necessary for simple scaling but for exibision purpose)
1315 : // 3. fit baseline to each spectrum and subtract it
1316 : // 4. put single spectrum (or a block of spectra) to out_column
1317 :
1318 836 : Block<Int> columns(1);
1319 418 : columns[0] = MS::DATA_DESC_ID;
1320 418 : prepare_for_process(in_column_name, out_ms_name, columns, false);
1321 418 : vi::VisibilityIterator2 *vi = sdh_->getVisIter();
1322 418 : vi->setRowBlocking(kNRowBlocking);
1323 418 : vi::VisBuffer2 *vb = vi->getVisBuffer();
1324 :
1325 418 : split_bloutputname(out_bloutput_name);
1326 418 : bool write_baseline_csv = (bloutputname_csv != "");
1327 418 : bool write_baseline_text = (bloutputname_text != "");
1328 418 : bool write_baseline_table = (bloutputname_table != "");
1329 836 : ofstream ofs_csv;
1330 836 : ofstream ofs_txt;
1331 418 : BaselineTable *bt = 0;
1332 :
1333 418 : if (write_baseline_csv) {
1334 195 : ofs_csv.open(bloutputname_csv.c_str());
1335 : }
1336 418 : if (write_baseline_text) {
1337 88 : ofs_txt.open(bloutputname_text.c_str(), std::ios_base::out | std::ios_base::app);
1338 : }
1339 418 : if (write_baseline_table) {
1340 53 : bt = new BaselineTable(vi->ms());
1341 : }
1342 :
1343 836 : Vector<Int> recspw;
1344 836 : Matrix<Int> recchan;
1345 836 : Vector<size_t> nchan;
1346 836 : Vector<Vector<Bool> > in_mask;
1347 836 : Vector<bool> nchan_set;
1348 418 : parse_spw(in_spw, recspw, recchan, nchan, in_mask, nchan_set);
1349 :
1350 836 : Vector<size_t> ctx_indices(nchan.nelements(), 0ul);
1351 836 : std::vector<int> blparam_eff_base;
1352 199072 : auto wn_ulimit_by_rejwn = [&](){
1353 198654 : return ((blparam_exclude.size() == 2) &&
1354 198654 : (blparam_exclude[1] == SinusoidWaveNumber_kUpperLimit)); };
1355 :
1356 836 : std::vector<float> weight(WeightIndex_kNum);
1357 418 : size_t const var_index = (sigma_value == "stddev") ? WeightIndex_kStddev : WeightIndex_kRms;
1358 :
1359 1261 : for (vi->originChunks(); vi->moreChunks(); vi->nextChunk()) {
1360 2462 : for (vi->origin(); vi->more(); vi->next()) {
1361 3238 : Vector<Int> scans = vb->scan();
1362 3238 : Vector<Double> times = vb->time();
1363 3238 : Vector<Int> beams = vb->feed1();
1364 3238 : Vector<Int> antennas = vb->antenna1();
1365 3238 : Vector<Int> data_spw = vb->spectralWindows();
1366 1619 : size_t const num_chan = static_cast<size_t>(vb->nChannels());
1367 1619 : size_t const num_pol = static_cast<size_t>(vb->nCorrelations());
1368 1619 : size_t const num_row = static_cast<size_t>(vb->nRows());
1369 3238 : Cube<Float> data_chunk(num_pol, num_chan, num_row, Array<Float>::uninitialized);
1370 3238 : SakuraAlignedArray<float> spec(num_chan);
1371 3238 : Cube<Bool> flag_chunk(num_pol, num_chan, num_row, Array<Bool>::uninitialized);
1372 3238 : SakuraAlignedArray<bool> flag(num_chan);
1373 3238 : SakuraAlignedArray<bool> mask(num_chan);
1374 3238 : SakuraAlignedArray<bool> mask_after_clipping(num_chan);
1375 1619 : float *spec_data = spec.data();
1376 1619 : bool *flag_data = flag.data();
1377 1619 : bool *mask_data = mask.data();
1378 1619 : bool *mask_after_clipping_data = mask_after_clipping.data();
1379 3238 : Matrix<Float> weight_matrix(num_pol, num_row, Array<Float>::uninitialized);
1380 :
1381 200339 : auto get_wavenumber_upperlimit = [&](){ return static_cast<int>(num_chan) / 2 - 1; };
1382 :
1383 1619 : uInt final_mask[num_pol];
1384 1619 : uInt final_mask_after_clipping[num_pol];
1385 1619 : final_mask[0] = 0;
1386 1619 : final_mask[1] = 0;
1387 1619 : final_mask_after_clipping[0] = 0;
1388 1619 : final_mask_after_clipping[1] = 0;
1389 :
1390 1619 : bool new_nchan = false;
1391 1619 : get_nchan_and_mask(recspw, data_spw, recchan, num_chan, nchan, in_mask, nchan_set, new_nchan);
1392 1619 : if (new_nchan) {
1393 418 : int blparam_max = blparam[blparam.size() - 1];
1394 418 : if (bltype == BaselineType_kSinusoid) {
1395 110 : int nwave_ulimit = get_wavenumber_upperlimit();
1396 110 : get_effective_nwave(blparam, blparam_exclude, nwave_ulimit, blparam_eff_base);
1397 104 : blparam_max = blparam_eff_base[blparam_eff_base.size() - 1];
1398 : }
1399 412 : if ((bltype != BaselineType_kSinusoid) || (!applyfft) || wn_ulimit_by_rejwn()) {
1400 331 : get_baseline_context(bltype,
1401 : static_cast<uint16_t>(blparam_max),
1402 : num_chan, nchan, nchan_set,
1403 : ctx_indices, bl_contexts);
1404 : }
1405 : } else {
1406 1201 : int last_nchan_set_idx = nchan_set.size() - 1;
1407 1549 : for (int i = nchan_set.size()-1; i >= 0; --i) {
1408 1549 : if (nchan_set[i]) break;
1409 348 : --last_nchan_set_idx;
1410 : }
1411 1201 : if (0 < last_nchan_set_idx) {
1412 431 : for (int i = 0; i < last_nchan_set_idx; ++i) {
1413 431 : if (nchan[i] == nchan[last_nchan_set_idx]) {
1414 431 : ctx_indices[last_nchan_set_idx] = ctx_indices[i];
1415 431 : break;
1416 : }
1417 : }
1418 : }
1419 : }
1420 :
1421 : // get data/flag cubes (npol*nchan*nrow) from VisBuffer
1422 1613 : get_data_cube_float(*vb, data_chunk);
1423 1613 : get_flag_cube(*vb, flag_chunk);
1424 :
1425 : // get weight matrix (npol*nrow) from VisBuffer
1426 1613 : if (update_weight) {
1427 23 : get_weight_matrix(*vb, weight_matrix);
1428 : }
1429 :
1430 : // loop over MS rows
1431 801260 : for (size_t irow = 0; irow < num_row; ++irow) {
1432 799647 : size_t idx = 0;
1433 800520 : for (size_t ispw = 0; ispw < recspw.nelements(); ++ispw) {
1434 800520 : if (data_spw[irow] == recspw[ispw]) {
1435 799647 : idx = ispw;
1436 799647 : break;
1437 : }
1438 : }
1439 :
1440 : //prepare variables for writing baseline table
1441 799647 : Array<Bool> apply_mtx(IPosition(2, num_pol, 1), true);
1442 799647 : Array<uInt> bltype_mtx(IPosition(2, num_pol, 1), (uInt)bltype);
1443 : //Array<Int> fpar_mtx(IPosition(2, num_pol, 1), (Int)blparam[blparam.size()-1]);
1444 799647 : std::vector<std::vector<size_t> > fpar_mtx_tmp(num_pol);
1445 799647 : std::vector<std::vector<double> > ffpar_mtx_tmp(num_pol);
1446 799647 : std::vector<std::vector<uInt> > masklist_mtx_tmp(num_pol);
1447 799647 : std::vector<std::vector<double> > coeff_mtx_tmp(num_pol);
1448 :
1449 799647 : Array<Float> rms_mtx(IPosition(2, num_pol, 1), (Float)0);
1450 799647 : Array<Float> cthres_mtx(IPosition(2, num_pol, 1), Array<Float>::uninitialized);
1451 799647 : Array<uInt> citer_mtx(IPosition(2, num_pol, 1), Array<uInt>::uninitialized);
1452 799647 : Array<Bool> uself_mtx(IPosition(2, num_pol, 1), Array<Bool>::uninitialized);
1453 799647 : Array<Float> lfthres_mtx(IPosition(2, num_pol, 1), Array<Float>::uninitialized);
1454 799647 : Array<uInt> lfavg_mtx(IPosition(2, num_pol, 1), Array<uInt>::uninitialized);
1455 799647 : Array<uInt> lfedge_mtx(IPosition(2, num_pol, 2), Array<uInt>::uninitialized);
1456 :
1457 799647 : size_t num_apply_true = 0;
1458 799647 : size_t num_fpar_max = 0;
1459 799647 : size_t num_ffpar_max = 0;
1460 799647 : size_t num_masklist_max = 0;
1461 799647 : size_t num_coeff_max = 0;
1462 :
1463 : // loop over polarization
1464 1602263 : for (size_t ipol = 0; ipol < num_pol; ++ipol) {
1465 : // get a channel mask from data cube
1466 : // (note that the variable 'mask' is flag in the next line
1467 : // actually, then it will be converted to real mask when
1468 : // taking AND with user-given mask info. this is just for
1469 : // saving memory usage...)
1470 802616 : get_flag_from_cube(flag_chunk, irow, ipol, num_chan, flag_data);
1471 : // skip spectrum if all channels flagged
1472 802616 : if (allchannels_flagged(num_chan, flag_data)) {
1473 3540 : apply_mtx[0][ipol] = false;
1474 3540 : continue;
1475 : }
1476 :
1477 : // convert flag to mask by taking logical NOT of flag
1478 : // and then operate logical AND with in_mask
1479 71645604 : for (size_t ichan = 0; ichan < num_chan; ++ichan) {
1480 70846528 : mask_data[ichan] = in_mask[idx][ichan] && (!(flag_data[ichan]));
1481 : }
1482 : // get a spectrum from data cube
1483 799076 : get_spectrum_from_cube(data_chunk, irow, ipol, num_chan, spec_data);
1484 : // line finding. get baseline mask (invert=true)
1485 799076 : if (linefinding) {
1486 3484 : findLineAndGetMask(num_chan, spec_data, mask_data, threshold,
1487 : avg_limit, minwidth, edge, true, mask_data);
1488 : }
1489 :
1490 799076 : std::vector<size_t> blparam_eff;
1491 :
1492 : size_t num_coeff;
1493 799076 : if (bltype == BaselineType_kSinusoid) {
1494 198610 : int nwave_ulimit = get_wavenumber_upperlimit();
1495 198610 : finalise_effective_nwave(blparam_eff_base, blparam_exclude,
1496 : nwave_ulimit,
1497 : num_chan, spec_data, mask_data,
1498 : applyfft, fftmethod, fftthresh,
1499 : blparam_eff);
1500 198610 : size_t blparam_eff_size = blparam_eff.size();
1501 198610 : if (blparam_eff[0] == 0) {
1502 198580 : num_coeff = blparam_eff_size * 2 - 1;
1503 : } else {
1504 30 : num_coeff = blparam_eff_size * 2;
1505 : }
1506 600466 : } else if (bltype == BaselineType_kCubicSpline) {
1507 198830 : blparam_eff.resize(1);
1508 198830 : blparam_eff[0] = blparam[blparam.size() - 1];
1509 198830 : num_coeff = blparam_eff[0] * 4;
1510 : } else { // poly, chebyshev
1511 401636 : blparam_eff.resize(1);
1512 401636 : blparam_eff[0] = blparam[blparam.size() - 1];
1513 401636 : status =
1514 401636 : LIBSAKURA_SYMBOL(GetNumberOfCoefficientsFloat)(bl_contexts[ctx_indices[idx]],
1515 : blparam_eff[0],
1516 : &num_coeff);
1517 401636 : check_sakura_status("sakura_GetNumberOfCoefficients", status);
1518 : }
1519 :
1520 : // Final check of the valid number of channels
1521 799076 : size_t num_min =
1522 198830 : (bltype == BaselineType_kCubicSpline) ? blparam[blparam.size()-1] + 3 : num_coeff;
1523 799076 : if (NValidMask(num_chan, mask_data) < num_min) {
1524 0 : flag_spectrum_in_cube(flag_chunk, irow, ipol);
1525 0 : apply_mtx[0][ipol] = false;
1526 : os << LogIO::WARN
1527 : << "Too few valid channels to fit. Skipping Antenna "
1528 0 : << antennas[irow] << ", Beam " << beams[irow] << ", SPW "
1529 0 : << data_spw[irow] << ", Pol " << ipol << ", Time "
1530 0 : << MVTime(times[irow] / 24. / 3600.).string(MVTime::YMD, 8)
1531 0 : << LogIO::POST;
1532 0 : continue;
1533 : }
1534 : // actual execution of single spectrum
1535 : float rms;
1536 799076 : if (write_baseline_text || write_baseline_csv || write_baseline_table) {
1537 798476 : num_apply_true++;
1538 :
1539 798476 : if (num_coeff_max < num_coeff) {
1540 795810 : num_coeff_max = num_coeff;
1541 : }
1542 1596952 : SakuraAlignedArray<double> coeff(num_coeff);
1543 798476 : double *coeff_data = coeff.data();
1544 :
1545 : //---GetBestFitBaselineCoefficientsFloat()...
1546 : //func0(ctx_indices[idx], num_chan, blparam_eff, spec_data, mask_data, num_coeff, coeff_data, mask_after_clipping_data, &rms);
1547 798476 : LIBSAKURA_SYMBOL(LSQFitContextFloat) *context = nullptr;
1548 798476 : if ((bltype != BaselineType_kSinusoid) || (!applyfft) || wn_ulimit_by_rejwn()) {
1549 600080 : context = bl_contexts[ctx_indices[idx]];
1550 : }
1551 798476 : func0(context, num_chan, blparam_eff, spec_data, mask_data, num_coeff, coeff_data, mask_after_clipping_data, &rms);
1552 :
1553 66729804 : for (size_t i = 0; i < num_chan; ++i) {
1554 65931328 : if (mask_data[i] == false) {
1555 10847408 : final_mask[ipol] += 1;
1556 : }
1557 65931328 : if (mask_after_clipping_data[i] == false) {
1558 10857281 : final_mask_after_clipping[ipol] += 1;
1559 : }
1560 : }
1561 :
1562 : //set_array_for_bltable(fpar_mtx_tmp)
1563 798476 : size_t num_fpar = blparam_eff.size();
1564 798476 : fpar_mtx_tmp[ipol].resize(num_fpar);
1565 798476 : if (num_fpar_max < num_fpar) {
1566 795810 : num_fpar_max = num_fpar;
1567 : }
1568 798476 : fpar_mtx_tmp[ipol].resize(num_fpar);
1569 1602420 : for (size_t ifpar = 0; ifpar < num_fpar; ++ifpar) {
1570 803944 : fpar_mtx_tmp[ipol][ifpar] = blparam_eff[ifpar];
1571 : }
1572 :
1573 : //---set_array_for_bltable(ffpar_mtx_tmp)
1574 798476 : func1(ipol, ffpar_mtx_tmp, num_ffpar_max);
1575 :
1576 : //set_array_for_bltable<double, Float>(ipol, num_coeff, coeff_data, coeff_mtx);
1577 798476 : coeff_mtx_tmp[ipol].resize(num_coeff);
1578 4222960 : for (size_t icoeff = 0; icoeff < num_coeff; ++icoeff) {
1579 3424484 : coeff_mtx_tmp[ipol][icoeff] = coeff_data[icoeff];
1580 : }
1581 :
1582 1596952 : Vector<uInt> masklist;
1583 798476 : get_masklist_from_mask(num_chan, mask_after_clipping_data, masklist);
1584 798476 : if (masklist.size() > num_masklist_max) {
1585 795827 : num_masklist_max = masklist.size();
1586 : }
1587 798476 : masklist_mtx_tmp[ipol].clear();
1588 3202316 : for (size_t imask = 0; imask < masklist.size(); ++imask) {
1589 2403840 : masklist_mtx_tmp[ipol].push_back(masklist[imask]);
1590 : }
1591 :
1592 : //---SubtractBaselineUsingCoefficientsFloat()...
1593 : //func2(ctx_indices[idx], num_chan, fpar_mtx_tmp[ipol], spec_data, num_coeff, coeff_data);
1594 798476 : func2(context, num_chan, fpar_mtx_tmp[ipol], spec_data, num_coeff, coeff_data);
1595 :
1596 798476 : rms_mtx[0][ipol] = rms;
1597 :
1598 798476 : cthres_mtx[0][ipol] = clip_threshold_sigma;
1599 798476 : citer_mtx[0][ipol] = (uInt)num_fitting_max - 1;
1600 798476 : uself_mtx[0][ipol] = false;
1601 798476 : lfthres_mtx[0][ipol] = 0.0;
1602 798476 : lfavg_mtx[0][ipol] = 0;
1603 2395428 : for (size_t iedge = 0; iedge < 2; ++iedge) {
1604 1596952 : lfedge_mtx[iedge][ipol] = 0;
1605 798476 : }
1606 : } else {
1607 : //---SubtractBaselineFloat()...
1608 : //func3(ctx_indices[idx], num_chan, blparam_eff, num_coeff, spec_data, mask_data, &rms);
1609 600 : LIBSAKURA_SYMBOL(LSQFitContextFloat) *context = nullptr;
1610 600 : if ((bltype != BaselineType_kSinusoid) || (!applyfft) || wn_ulimit_by_rejwn()) {
1611 432 : context = bl_contexts[ctx_indices[idx]];
1612 : }
1613 600 : func3(context, num_chan, blparam_eff, num_coeff, spec_data, mask_data, mask_after_clipping_data, &rms);
1614 : }
1615 : // set back a spectrum to data cube
1616 799076 : if (do_subtract) {
1617 799016 : set_spectrum_to_cube(data_chunk, irow, ipol, num_chan, spec_data);
1618 : }
1619 :
1620 799076 : if (update_weight) {
1621 4680 : compute_weight(num_chan, spec_data, mask_after_clipping_data, weight);
1622 4680 : set_weight_to_matrix(weight_matrix, irow, ipol, weight.at(var_index));
1623 : }
1624 :
1625 : } // end of polarization loop
1626 :
1627 : // output results of fitting
1628 799647 : if (num_apply_true == 0) continue;
1629 :
1630 1591620 : Array<Int> fpar_mtx(IPosition(2, num_pol, num_fpar_max),
1631 : Array<Int>::uninitialized);
1632 795810 : set_matrix_for_bltable<size_t, Int>(num_pol, num_fpar_max,
1633 : fpar_mtx_tmp, fpar_mtx);
1634 1591620 : Array<Float> ffpar_mtx(IPosition(2, num_pol, num_ffpar_max),
1635 : Array<Float>::uninitialized);
1636 795810 : set_matrix_for_bltable<double, Float>(num_pol, num_ffpar_max,
1637 : ffpar_mtx_tmp, ffpar_mtx);
1638 1591620 : Array<uInt> masklist_mtx(IPosition(2, num_pol, num_masklist_max),
1639 : Array<uInt>::uninitialized);
1640 795810 : set_matrix_for_bltable<uInt, uInt>(num_pol, num_masklist_max,
1641 : masklist_mtx_tmp, masklist_mtx);
1642 1591620 : Array<Float> coeff_mtx(IPosition(2, num_pol, num_coeff_max),
1643 : Array<Float>::uninitialized);
1644 795810 : set_matrix_for_bltable<double, Float>(num_pol, num_coeff_max,
1645 : coeff_mtx_tmp, coeff_mtx);
1646 1591620 : Matrix<uInt> masklist_mtx2 = masklist_mtx;
1647 1591620 : Matrix<Bool> apply_mtx2 = apply_mtx;
1648 :
1649 795810 : if (write_baseline_table) {
1650 378 : bt->appenddata((uInt)scans[irow], (uInt)beams[irow], (uInt)antennas[irow],
1651 189 : (uInt)data_spw[irow], 0, times[irow],
1652 : apply_mtx, bltype_mtx, fpar_mtx, ffpar_mtx, masklist_mtx,
1653 : coeff_mtx, rms_mtx, (uInt)num_chan, cthres_mtx, citer_mtx,
1654 : uself_mtx, lfthres_mtx, lfavg_mtx, lfedge_mtx);
1655 : }
1656 :
1657 795810 : if (write_baseline_text) {
1658 7349 : for (size_t ipol = 0; ipol < num_pol; ++ipol) {
1659 4896 : if (apply_mtx2(ipol, 0) == false) continue;
1660 :
1661 4896 : ofs_txt << "Scan" << '[' << (uInt)scans[irow] << ']' << ' '
1662 4896 : << "Beam" << '[' << (uInt)beams[irow] << ']' << ' '
1663 4896 : << "Spw" << '[' << (uInt)data_spw[irow] << ']' << ' '
1664 4896 : << "Pol" << '[' << ipol << ']' << ' '
1665 9792 : << "Time" << '[' << MVTime(times[irow]/ 24. / 3600.).string(MVTime::YMD, 8) << ']'
1666 4896 : << endl;
1667 4896 : ofs_txt << endl;
1668 4896 : ofs_txt << "Fitter range = " << '[';
1669 :
1670 9865 : for (size_t imasklist = 0; imasklist < num_masklist_max/2; ++imasklist) {
1671 4969 : if (imasklist == 0) {
1672 4896 : ofs_txt << '[' << masklist_mtx2(ipol, 2 * imasklist) << ';'
1673 4896 : << masklist_mtx2(ipol, 2 * imasklist + 1) << ']';
1674 : }
1675 4969 : if (imasklist >= 1
1676 5041 : && (0 != masklist_mtx2(ipol, 2 * imasklist)
1677 72 : && 0 != masklist_mtx2(ipol, 2 * imasklist + 1))) {
1678 72 : ofs_txt << ",[" << masklist_mtx2(ipol, 2 * imasklist) << ','
1679 72 : << masklist_mtx2(ipol, 2 * imasklist + 1) << ']';
1680 : }
1681 : }
1682 :
1683 4896 : ofs_txt << ']' << endl;
1684 4896 : ofs_txt << endl;
1685 :
1686 14688 : Matrix<uInt> bltype_mtx2 = bltype_mtx[0][ipol];
1687 14688 : Matrix<Int> fpar_mtx2 = fpar_mtx[0][ipol];
1688 14688 : Matrix<Float> rms_mtx2 = rms_mtx[0][ipol];
1689 9792 : string bltype_name;
1690 :
1691 9792 : string blparam_name ="order";
1692 4896 : if (bltype_mtx2(0, 0) == (uInt)0) {
1693 4811 : bltype_name = "poly";
1694 85 : } else if (bltype_mtx2(0, 0) == (uInt)1) {
1695 1 : bltype_name = "chebyshev";
1696 84 : } else if (bltype_mtx2(0, 0) == (uInt)2) {
1697 2 : blparam_name = "npiece";
1698 2 : bltype_name = "cspline";
1699 82 : } else if (bltype_mtx2(0, 0) == (uInt)3) {
1700 82 : blparam_name = "nwave";
1701 82 : bltype_name = "sinusoid";
1702 : }
1703 :
1704 : ofs_txt << "Baseline parameters Function = "
1705 4896 : << bltype_name << " " << blparam_name << " = ";
1706 9792 : Matrix<Int> fpar_mtx3 = fpar_mtx;
1707 4896 : if (bltype_mtx2(0,0) == (uInt)3) {
1708 504 : for (size_t num = 0; num < num_fpar_max; ++num) {
1709 422 : ofs_txt << fpar_mtx3(ipol, num) << " ";
1710 : }
1711 82 : ofs_txt << endl;
1712 : } else {
1713 4814 : ofs_txt << fpar_mtx2(0, 0) << endl;
1714 : }
1715 :
1716 4896 : ofs_txt << endl;
1717 4896 : ofs_txt << "Results of baseline fit" << endl;
1718 4896 : ofs_txt << endl;
1719 9792 : Matrix<Float> coeff_mtx2 = coeff_mtx;
1720 :
1721 4896 : if (bltype_mtx2(0,0) == (uInt)0 || bltype_mtx2(0,0) == (uInt)1 || bltype_mtx2(0,0) == (uInt)2){
1722 33518 : for (size_t icoeff = 0; icoeff < num_coeff_max; ++icoeff) {
1723 28704 : ofs_txt << "p" << icoeff << " = " << setprecision(8) << coeff_mtx2(ipol, icoeff) << " ";
1724 : }
1725 82 : } else if (bltype_mtx2(0,0) == (uInt)3) {
1726 82 : size_t wn=0;
1727 164 : string c_s ="s";
1728 : //if (blparam[0] == 0) {
1729 82 : if (fpar_mtx3(ipol, wn) == 0) {
1730 52 : ofs_txt << "c" << fpar_mtx3(ipol, wn) << " = " <<setw(13)<<left<< setprecision(8) << coeff_mtx2(ipol, 0) << " ";
1731 52 : wn = 1;
1732 : //for (size_t icoeff = 1; icoeff < num_coeff_max; ++icoeff) {
1733 232 : for (size_t icoeff = 1; icoeff < coeff_mtx_tmp[ipol].size(); ++icoeff) {
1734 180 : ofs_txt << c_s << fpar_mtx3(ipol, wn) << " = " <<setw(13)<<left<< setprecision(8) << coeff_mtx2(ipol, icoeff) << " ";
1735 180 : c_s == "s" ? (c_s = "c") : (c_s = "s");
1736 180 : if (icoeff % 2 == 0) {
1737 90 : ++wn;
1738 : }
1739 : }
1740 : } else {
1741 30 : wn = 0;
1742 : //for (size_t icoeff = 0; icoeff < num_coeff_max; ++icoeff) {
1743 590 : for (size_t icoeff = 0; icoeff < coeff_mtx_tmp[ipol].size(); ++icoeff) {
1744 560 : ofs_txt << c_s << fpar_mtx3(ipol, wn) << " = " <<setw(13)<<left<< setprecision(8) << coeff_mtx2(ipol, icoeff) << " ";
1745 560 : c_s == "s" ? (c_s = "c") : (c_s = "s");
1746 560 : if (icoeff % 2 != 0) {
1747 280 : ++wn;
1748 : }
1749 : }
1750 : }
1751 : }
1752 :
1753 4896 : ofs_txt << endl;
1754 4896 : ofs_txt << endl;
1755 4896 : ofs_txt << "rms = ";
1756 4896 : ofs_txt << setprecision(8) << rms_mtx2(0, 0) << endl;
1757 4896 : ofs_txt << endl;
1758 4896 : ofs_txt << "Number of clipped channels = "
1759 4896 : << final_mask_after_clipping[ipol] - final_mask[ipol] << endl;
1760 4896 : ofs_txt << endl;
1761 4896 : ofs_txt << "------------------------------------------------------"
1762 4896 : << endl;
1763 4896 : ofs_txt << endl;
1764 : }
1765 : }
1766 :
1767 795810 : if (write_baseline_csv) {
1768 1586475 : for (size_t ipol = 0; ipol < num_pol; ++ipol) {
1769 793278 : if (apply_mtx2(ipol, 0) == false) continue;
1770 :
1771 793278 : ofs_csv << (uInt)scans[irow] << ',' << (uInt)beams[irow] << ','
1772 793278 : << (uInt)data_spw[irow] << ',' << ipol << ','
1773 793278 : << setprecision(12) << times[irow] << ',';
1774 793278 : ofs_csv << '[';
1775 1989922 : for (size_t imasklist = 0; imasklist < num_masklist_max / 2; ++imasklist) {
1776 1196644 : if (imasklist == 0) {
1777 793278 : ofs_csv << '[' << masklist_mtx2(ipol, 2 * imasklist) << ';'
1778 793278 : << masklist_mtx2(ipol, 2 * imasklist + 1) << ']';
1779 : }
1780 1196644 : if (imasklist >= 1
1781 1599981 : && (0 != masklist_mtx2(ipol, 2 * imasklist)
1782 403337 : && 0 != masklist_mtx2(ipol, 2 * imasklist + 1))) {
1783 403337 : ofs_csv << ";[" << masklist_mtx2(ipol, 2 * imasklist) << ';'
1784 403337 : << masklist_mtx2(ipol, 2 * imasklist + 1) << ']';
1785 : }
1786 : }
1787 793278 : ofs_csv << ']' << ',';
1788 2379834 : Matrix<uInt> bltype_mtx2 = bltype_mtx[0][ipol];
1789 1586556 : string bltype_name;
1790 793278 : if (bltype_mtx2(0, 0) == (uInt)0) {
1791 198268 : bltype_name = "poly";
1792 595010 : } else if (bltype_mtx2(0, 0) == (uInt)1) {
1793 198148 : bltype_name = "chebyshev";
1794 396862 : } else if (bltype_mtx2(0, 0) == (uInt)2) {
1795 198556 : bltype_name = "cspline";
1796 198306 : } else if (bltype_mtx2(0, 0) == (uInt)3) {
1797 198306 : bltype_name = "sinusoid";
1798 : }
1799 :
1800 1586556 : Matrix<Int> fpar_mtx2 = fpar_mtx;
1801 793278 : if (bltype_mtx2(0, 0) == (uInt)3) {
1802 198306 : ofs_csv << bltype_name.c_str() << ',' << fpar_mtx2(ipol, 0);
1803 203644 : for (size_t i = 1; i < num_fpar_max; ++i) {
1804 5338 : ofs_csv << ';' << fpar_mtx2(ipol, i);
1805 : }
1806 198306 : ofs_csv << ',';
1807 : } else {
1808 594972 : ofs_csv << bltype_name.c_str() << ',' << fpar_mtx2(ipol, 0)
1809 594972 : << ',';
1810 : }
1811 :
1812 1586556 : Matrix<Float> coeff_mtx2 = coeff_mtx;
1813 793278 : if (bltype_mtx2(0, 0) == (uInt)3) {
1814 407290 : for (size_t icoeff = 0; icoeff < coeff_mtx_tmp[ipol].size(); ++icoeff) {
1815 208984 : ofs_csv << setprecision(8) << coeff_mtx2(ipol, icoeff) << ',';
1816 : }
1817 : } else {
1818 3779476 : for (size_t icoeff = 0; icoeff < num_coeff_max; ++icoeff) {
1819 3184504 : ofs_csv << setprecision(8) << coeff_mtx2(ipol, icoeff) << ',';
1820 : }
1821 : }
1822 :
1823 1586556 : Matrix<Float> rms_mtx2 = rms_mtx;
1824 793278 : ofs_csv << setprecision(8) << rms_mtx2(ipol, 0) << ',';
1825 793278 : ofs_csv << final_mask_after_clipping[ipol] - final_mask[ipol];
1826 793278 : ofs_csv << endl;
1827 : }
1828 : }
1829 : } // end of chunk row loop
1830 : // write back data cube to VisBuffer
1831 1613 : if (do_subtract) {
1832 1583 : if (update_weight) {
1833 23 : sdh_->fillCubeToOutputMs(vb, data_chunk, &flag_chunk, &weight_matrix);
1834 : } else {
1835 1560 : sdh_->fillCubeToOutputMs(vb, data_chunk, &flag_chunk);
1836 : }
1837 : }
1838 : } // end of vi loop
1839 : } // end of chunk loop
1840 :
1841 412 : if (write_baseline_table) {
1842 53 : bt->save(bloutputname_table);
1843 53 : delete bt;
1844 : }
1845 412 : if (write_baseline_csv) {
1846 195 : ofs_csv.close();
1847 : }
1848 412 : if (write_baseline_text) {
1849 82 : ofs_txt.close();
1850 : }
1851 :
1852 412 : finalize_process();
1853 412 : destroy_baseline_contexts(bl_contexts);
1854 :
1855 : //double tend = gettimeofday_sec();
1856 : //std::cout << "Elapsed time = " << (tend - tstart) << " sec." << std::endl;
1857 412 : }
1858 :
1859 : ////////////////////////////////////////////////////////////////////////
1860 : ///// Atcual processing functions
1861 : ////////////////////////////////////////////////////////////////////////
1862 :
1863 : //Subtract baseline using normal or Chebyshev polynomials
1864 226 : void SingleDishMS::subtractBaseline(string const& in_column_name,
1865 : string const& out_ms_name,
1866 : string const& out_bloutput_name,
1867 : bool const& do_subtract,
1868 : string const& in_spw,
1869 : bool const& update_weight,
1870 : string const& sigma_value,
1871 : string const& blfunc,
1872 : int const order,
1873 : float const clip_threshold_sigma,
1874 : int const num_fitting_max,
1875 : bool const linefinding,
1876 : float const threshold,
1877 : int const avg_limit,
1878 : int const minwidth,
1879 : vector<int> const& edge) {
1880 452 : vector<int> order_vect;
1881 226 : order_vect.push_back(order);
1882 452 : vector<int> blparam_exclude_dummy;
1883 226 : blparam_exclude_dummy.clear();
1884 :
1885 678 : LogIO os(_ORIGIN);
1886 : os << "Fitting and subtracting polynomial baseline order = " << order
1887 226 : << LogIO::POST;
1888 226 : if (order < 0) {
1889 0 : throw(AipsError("order must be positive or zero."));
1890 : }
1891 :
1892 : LIBSAKURA_SYMBOL(Status) status;
1893 : LIBSAKURA_SYMBOL(LSQFitStatus) bl_status;
1894 452 : std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> bl_contexts;
1895 226 : bl_contexts.clear();
1896 226 : size_t bltype = BaselineType_kPolynomial;
1897 226 : string blfunc_lower = blfunc;
1898 : std::transform(
1899 : blfunc_lower.begin(),
1900 : blfunc_lower.end(),
1901 : blfunc_lower.begin(),
1902 1134 : [](unsigned char c) {return std::tolower(c);}
1903 226 : );
1904 226 : if (blfunc_lower == "chebyshev") {
1905 46 : bltype = BaselineType_kChebyshev;
1906 : }
1907 :
1908 452 : doSubtractBaseline(in_column_name,
1909 : out_ms_name,
1910 : out_bloutput_name,
1911 : do_subtract,
1912 : in_spw,
1913 : update_weight,
1914 : sigma_value,
1915 : status,
1916 : bl_contexts,
1917 : bltype,
1918 : order_vect,
1919 : blparam_exclude_dummy,
1920 226 : false,
1921 : "",
1922 : "",
1923 : clip_threshold_sigma,
1924 : num_fitting_max,
1925 : linefinding,
1926 : threshold,
1927 : avg_limit,
1928 : minwidth,
1929 : edge,
1930 : [&](LIBSAKURA_SYMBOL(LSQFitContextFloat) const *context,
1931 : size_t const num_chan, std::vector<size_t> const &/*nwave*/,
1932 : float *spec, bool const *mask, size_t const /*num_coeff*/, double *coeff,
1933 401372 : bool *mask_after_clipping, float *rms){
1934 401372 : status = LIBSAKURA_SYMBOL(LSQFitPolynomialFloat)(
1935 802744 : context, static_cast<uint16_t>(order_vect[0]),
1936 802744 : num_chan, spec, mask, clip_threshold_sigma, num_fitting_max,
1937 401372 : order_vect[0] + 1, coeff, nullptr, nullptr, mask_after_clipping, rms, &bl_status);
1938 401372 : check_sakura_status("sakura_LSQFitPolynomialFloat", status);
1939 401372 : if (bl_status != LIBSAKURA_SYMBOL(LSQFitStatus_kOK)) {
1940 0 : throw(AipsError("baseline fitting isn't successful."));
1941 : }
1942 401372 : },
1943 401372 : [&](size_t ipol, std::vector<std::vector<double> > &ffpar_mtx_tmp, size_t &num_ffpar_max) {
1944 401372 : size_t num_ffpar = get_num_coeff_bloutput(bltype, 0, num_ffpar_max);
1945 401372 : ffpar_mtx_tmp[ipol].resize(num_ffpar);
1946 401372 : },
1947 : [&](LIBSAKURA_SYMBOL(LSQFitContextFloat) const *context,
1948 : size_t const num_chan, std::vector<size_t> const &/*nwave*/,
1949 401372 : float *spec, size_t const /*num_coeff*/, double *coeff){
1950 1204116 : status = LIBSAKURA_SYMBOL(SubtractPolynomialFloat)(
1951 401372 : context, num_chan, spec, order_vect[0] + 1, coeff, spec);
1952 401372 : check_sakura_status("sakura_SubtractPolynomialFloat", status);},
1953 : [&](LIBSAKURA_SYMBOL(LSQFitContextFloat) const *context,
1954 : size_t const num_chan, std::vector<size_t> const &/*nwave*/,
1955 264 : size_t const /*num_coeff*/, float *spec, bool const *mask, bool *mask_after_clipping, float *rms){
1956 264 : status = LIBSAKURA_SYMBOL(LSQFitPolynomialFloat)(
1957 528 : context, static_cast<uint16_t>(order_vect[0]),
1958 528 : num_chan, spec, mask, clip_threshold_sigma, num_fitting_max,
1959 264 : order_vect[0] + 1, nullptr, nullptr, spec, mask_after_clipping, rms, &bl_status);
1960 264 : check_sakura_status("sakura_LSQFitPolynomialFloat", status);
1961 264 : if (bl_status != LIBSAKURA_SYMBOL(LSQFitStatus_kOK)) {
1962 0 : throw(AipsError("baseline fitting isn't successful."));
1963 : }
1964 264 : },
1965 : os
1966 : );
1967 226 : }
1968 :
1969 : //Subtract baseline using natural cubic spline
1970 82 : void SingleDishMS::subtractBaselineCspline(string const& in_column_name,
1971 : string const& out_ms_name,
1972 : string const& out_bloutput_name,
1973 : bool const& do_subtract,
1974 : string const& in_spw,
1975 : bool const& update_weight,
1976 : string const& sigma_value,
1977 : int const npiece,
1978 : float const clip_threshold_sigma,
1979 : int const num_fitting_max,
1980 : bool const linefinding,
1981 : float const threshold,
1982 : int const avg_limit,
1983 : int const minwidth,
1984 : vector<int> const& edge) {
1985 164 : vector<int> npiece_vect;
1986 82 : npiece_vect.push_back(npiece);
1987 164 : vector<int> blparam_exclude_dummy;
1988 82 : blparam_exclude_dummy.clear();
1989 :
1990 246 : LogIO os(_ORIGIN);
1991 : os << "Fitting and subtracting cubic spline baseline npiece = " << npiece
1992 82 : << LogIO::POST;
1993 82 : if (npiece <= 0) {
1994 0 : throw(AipsError("npiece must be positive."));
1995 : }
1996 :
1997 : LIBSAKURA_SYMBOL(Status) status;
1998 : LIBSAKURA_SYMBOL(LSQFitStatus) bl_status;
1999 164 : std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> bl_contexts;
2000 82 : bl_contexts.clear();
2001 82 : size_t const bltype = BaselineType_kCubicSpline;
2002 82 : SakuraAlignedArray<size_t> boundary(npiece+1);
2003 82 : size_t *boundary_data = boundary.data();
2004 :
2005 164 : doSubtractBaseline(in_column_name,
2006 : out_ms_name,
2007 : out_bloutput_name,
2008 : do_subtract,
2009 : in_spw,
2010 : update_weight,
2011 : sigma_value,
2012 : status,
2013 : bl_contexts,
2014 : bltype,
2015 : npiece_vect,
2016 : blparam_exclude_dummy,
2017 82 : false,
2018 : "",
2019 : "",
2020 : clip_threshold_sigma,
2021 : num_fitting_max,
2022 : linefinding,
2023 : threshold,
2024 : avg_limit,
2025 : minwidth,
2026 : edge,
2027 : [&](LIBSAKURA_SYMBOL(LSQFitContextFloat) const *context,
2028 : size_t const num_chan, std::vector<size_t> const &/*nwave*/,
2029 : float *spec, bool const *mask, size_t const /*num_coeff*/, double *coeff,
2030 198662 : bool *mask_after_clipping, float *rms) {
2031 198662 : status = LIBSAKURA_SYMBOL(LSQFitCubicSplineFloat)(
2032 198662 : context, static_cast<uint16_t>(npiece_vect[0]),
2033 397324 : num_chan, spec, mask, clip_threshold_sigma, num_fitting_max,
2034 : reinterpret_cast<double (*)[4]>(coeff), nullptr, nullptr,
2035 198662 : mask_after_clipping, rms, boundary_data, &bl_status);
2036 198662 : check_sakura_status("sakura_LSQFitCubicSplineFloat", status);
2037 198662 : if (bl_status != LIBSAKURA_SYMBOL(LSQFitStatus_kOK)) {
2038 0 : throw(AipsError("baseline fitting isn't successful."));
2039 : }
2040 198662 : },
2041 198662 : [&](size_t ipol, std::vector<std::vector<double> > &ffpar_mtx_tmp, size_t &num_ffpar_max) {
2042 198662 : size_t num_ffpar = get_num_coeff_bloutput(
2043 198662 : bltype, npiece_vect[0], num_ffpar_max);
2044 198662 : ffpar_mtx_tmp[ipol].resize(num_ffpar);
2045 992435 : for (size_t ipiece = 0; ipiece < num_ffpar; ++ipiece) {
2046 793773 : ffpar_mtx_tmp[ipol][ipiece] = boundary_data[ipiece];
2047 : }
2048 198662 : },
2049 : [&](LIBSAKURA_SYMBOL(LSQFitContextFloat) const *context,
2050 : size_t const num_chan, std::vector<size_t> const &/*nwave*/,
2051 198662 : float *spec, size_t const /*num_coeff*/, double *coeff) {
2052 198662 : status = LIBSAKURA_SYMBOL(SubtractCubicSplineFloat)(
2053 198662 : context, num_chan, spec, npiece_vect[0],
2054 198662 : reinterpret_cast<double (*)[4]>(coeff), boundary_data, spec);
2055 198662 : check_sakura_status("sakura_SubtractCubicSplineFloat", status);},
2056 : [&](LIBSAKURA_SYMBOL(LSQFitContextFloat) const *context,
2057 : size_t const num_chan, std::vector<size_t> const &/*nwave*/,
2058 168 : size_t const /*num_coeff*/, float *spec, bool const *mask, bool *mask_after_clipping, float *rms) {
2059 168 : status = LIBSAKURA_SYMBOL(LSQFitCubicSplineFloat)(
2060 168 : context, static_cast<uint16_t>(npiece_vect[0]),
2061 336 : num_chan, spec, mask, clip_threshold_sigma, num_fitting_max,
2062 168 : nullptr, nullptr, spec, mask_after_clipping, rms, boundary_data, &bl_status);
2063 168 : check_sakura_status("sakura_LSQFitCubicSplineFloat", status);
2064 168 : if (bl_status != LIBSAKURA_SYMBOL(LSQFitStatus_kOK)) {
2065 0 : throw(AipsError("baseline fitting isn't successful."));
2066 : }
2067 168 : },
2068 : os
2069 : );
2070 82 : }
2071 :
2072 :
2073 114 : void SingleDishMS::subtractBaselineSinusoid(string const& in_column_name,
2074 : string const& out_ms_name,
2075 : string const& out_bloutput_name,
2076 : bool const& do_subtract,
2077 : string const& in_spw,
2078 : bool const& update_weight,
2079 : string const& sigma_value,
2080 : string const& addwn0,
2081 : string const& rejwn0,
2082 : bool const applyfft,
2083 : string const fftmethod,
2084 : string const fftthresh,
2085 : float const clip_threshold_sigma,
2086 : int const num_fitting_max,
2087 : bool const linefinding,
2088 : float const threshold,
2089 : int const avg_limit,
2090 : int const minwidth,
2091 : vector<int> const& edge) {
2092 114 : char const delim = ',';
2093 228 : vector<int> addwn = string_to_list(addwn0, delim);
2094 228 : vector<int> rejwn = string_to_list(rejwn0, delim);
2095 :
2096 342 : LogIO os(_ORIGIN);
2097 114 : os << "Fitting and subtracting sinusoid baseline with wave numbers " << addwn0 << LogIO::POST;
2098 114 : if (addwn.size() == 0) {
2099 4 : throw(AipsError("addwn must contain at least one element."));
2100 : }
2101 :
2102 : LIBSAKURA_SYMBOL(Status) status;
2103 : LIBSAKURA_SYMBOL(LSQFitStatus) bl_status;
2104 116 : std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> bl_contexts;
2105 110 : bl_contexts.clear();
2106 110 : LIBSAKURA_SYMBOL(LSQFitContextFloat) *context = nullptr;
2107 110 : size_t bltype = BaselineType_kSinusoid;
2108 :
2109 595524 : auto wn_ulimit_by_rejwn = [&](){
2110 595548 : return ((rejwn.size() == 2) &&
2111 595548 : (rejwn[1] == SinusoidWaveNumber_kUpperLimit)); };
2112 595662 : auto par_spectrum_context = [&](){
2113 595662 : return (applyfft && !wn_ulimit_by_rejwn());
2114 110 : };
2115 : auto prepare_context = [&](LIBSAKURA_SYMBOL(LSQFitContextFloat) const *context0,
2116 198610 : size_t const num_chan, std::vector<size_t> const &nwave){
2117 198610 : if (par_spectrum_context()) {
2118 198564 : status = LIBSAKURA_SYMBOL(CreateLSQFitContextSinusoidFloat)(
2119 198564 : static_cast<uint16_t>(nwave[nwave.size()-1]),
2120 198564 : num_chan, &context);
2121 198564 : check_sakura_status("sakura_CreateLSQFitContextSinusoidFloat", status);
2122 : } else {
2123 46 : context = const_cast<LIBSAKURA_SYMBOL(LSQFitContextFloat) *>(context0);
2124 : }
2125 198610 : };
2126 198610 : auto clear_context = [&](){
2127 198610 : if (par_spectrum_context()) {
2128 198564 : status = LIBSAKURA_SYMBOL(DestroyLSQFitContextFloat)(context);
2129 198564 : check_sakura_status("sakura_DestoyBaselineContextFloat", status);
2130 198564 : context = nullptr;
2131 : }
2132 198610 : };
2133 :
2134 116 : doSubtractBaseline(in_column_name,
2135 : out_ms_name,
2136 : out_bloutput_name,
2137 : do_subtract,
2138 : in_spw,
2139 : update_weight,
2140 : sigma_value,
2141 : status,
2142 : bl_contexts,
2143 : bltype,
2144 : addwn,
2145 : rejwn,
2146 : applyfft,
2147 : fftmethod,
2148 : fftthresh,
2149 : clip_threshold_sigma,
2150 : num_fitting_max,
2151 : linefinding,
2152 : threshold,
2153 : avg_limit,
2154 : minwidth,
2155 : edge,
2156 : [&](LIBSAKURA_SYMBOL(LSQFitContextFloat) const *context0,
2157 : size_t const num_chan, std::vector<size_t> const &nwave,
2158 : float *spec, bool const *mask, size_t const num_coeff, double *coeff,
2159 198442 : bool *mask_after_clipping, float *rms) {
2160 198442 : prepare_context(context0, num_chan, nwave);
2161 198442 : status = LIBSAKURA_SYMBOL(LSQFitSinusoidFloat)(
2162 198442 : context, nwave.size(), &nwave[0],
2163 396884 : num_chan, spec, mask, clip_threshold_sigma, num_fitting_max,
2164 198442 : num_coeff, coeff, nullptr, nullptr, mask_after_clipping, rms, &bl_status);
2165 198442 : check_sakura_status("sakura_LSQFitSinusoidFloat", status);
2166 198442 : check_baseline_status(bl_status);
2167 198442 : },
2168 198442 : [&](size_t ipol, std::vector<std::vector<double> > &ffpar_mtx_tmp, size_t &num_ffpar_max) {
2169 198442 : size_t num_ffpar = get_num_coeff_bloutput(bltype, addwn.size(), num_ffpar_max);
2170 198442 : ffpar_mtx_tmp[ipol].resize(num_ffpar);
2171 198442 : },
2172 : [&](LIBSAKURA_SYMBOL(LSQFitContextFloat) const *context0,
2173 : size_t const num_chan, std::vector<size_t> const &nwave,
2174 198442 : float *spec, size_t num_coeff, double *coeff) {
2175 198442 : if (!par_spectrum_context()) {
2176 46 : context = const_cast<LIBSAKURA_SYMBOL(LSQFitContextFloat) *>(context0);
2177 : }
2178 198442 : status = LIBSAKURA_SYMBOL(SubtractSinusoidFloat)(
2179 198442 : context, num_chan, spec, nwave.size(), &nwave[0],
2180 : num_coeff, coeff, spec);
2181 198442 : check_sakura_status("sakura_SubtractSinusoidFloat", status);
2182 198442 : clear_context();
2183 198442 : },
2184 : [&](LIBSAKURA_SYMBOL(LSQFitContextFloat) const *context0,
2185 : size_t const num_chan, std::vector<size_t> const &nwave,
2186 168 : size_t const num_coeff, float *spec, bool const *mask, bool *mask_after_clipping, float *rms) {
2187 168 : prepare_context(context0, num_chan, nwave);
2188 168 : status = LIBSAKURA_SYMBOL(LSQFitSinusoidFloat)(
2189 168 : context, nwave.size(), &nwave[0],
2190 336 : num_chan, spec, mask, clip_threshold_sigma, num_fitting_max,
2191 168 : num_coeff, nullptr, nullptr, spec, mask_after_clipping, rms, &bl_status);
2192 168 : check_sakura_status("sakura_LSQFitSinusoidFloat", status);
2193 168 : check_baseline_status(bl_status);
2194 168 : clear_context();
2195 168 : },
2196 : os
2197 : );
2198 104 : }
2199 :
2200 : // Apply baseline table to MS
2201 10 : void SingleDishMS::applyBaselineTable(string const& in_column_name,
2202 : string const& in_bltable_name,
2203 : string const& in_spw,
2204 : bool const& update_weight,
2205 : string const& sigma_value,
2206 : string const& out_ms_name) {
2207 30 : LogIO os(_ORIGIN);
2208 10 : os << "Apply baseline table " << in_bltable_name << " to MS. " << LogIO::POST;
2209 :
2210 10 : if (in_bltable_name == "") {
2211 0 : throw(AipsError("baseline table is not given."));
2212 : }
2213 :
2214 : // parse fitting parameters in the text file
2215 20 : BLTableParser parser(in_bltable_name);
2216 20 : std::vector<size_t> baseline_types = parser.get_function_types();
2217 20 : map<size_t const, uint16_t> max_orders;
2218 30 : for (size_t i = 0; i < baseline_types.size(); ++i) {
2219 20 : max_orders[baseline_types[i]] = parser.get_max_order(baseline_types[i]);
2220 : }
2221 : { //DEBUG output
2222 10 : os << LogIO::DEBUG1 << "spw ID = " << in_spw << LogIO::POST;
2223 10 : os << LogIO::DEBUG1 << "Baseline Types = " << baseline_types << LogIO::POST;
2224 10 : os << LogIO::DEBUG1 << "Max Orders:" << LogIO::POST;
2225 10 : map<size_t const, uint16_t>::iterator iter = max_orders.begin();
2226 30 : while (iter != max_orders.end()) {
2227 40 : os << LogIO::DEBUG1 << "- type " << (*iter).first << ": "
2228 20 : << (*iter).second << LogIO::POST;
2229 20 : ++iter;
2230 : }
2231 : }
2232 :
2233 20 : Block<Int> columns(1);
2234 10 : columns[0] = MS::DATA_DESC_ID; // (CAS-9918, 2017/4/27 WK) //columns[0] = MS::TIME;
2235 10 : prepare_for_process(in_column_name, out_ms_name, columns, false);
2236 10 : vi::VisibilityIterator2 *vi = sdh_->getVisIter();
2237 10 : vi->setRowBlocking(kNRowBlocking);
2238 10 : vi::VisBuffer2 *vb = vi->getVisBuffer();
2239 :
2240 20 : Vector<Int> recspw;
2241 20 : Matrix<Int> recchan;
2242 20 : Vector<size_t> nchan;
2243 20 : Vector<Vector<Bool> > in_mask;
2244 20 : Vector<bool> nchan_set;
2245 10 : parse_spw(in_spw, recspw, recchan, nchan, in_mask, nchan_set);
2246 :
2247 : // Baseline Contexts reservoir
2248 : // key: BaselineType
2249 : // value: a vector of Sakura_BaselineContextFloat for various nchans
2250 20 : Vector<size_t> ctx_indices(nchan.nelements(), 0ul);
2251 20 : map<size_t const, std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> > context_reservoir;
2252 10 : map<size_t const, uint16_t>::iterator iter = max_orders.begin();
2253 30 : while (iter != max_orders.end()) {
2254 20 : context_reservoir[(*iter).first] = std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *>();
2255 20 : ++iter;
2256 : }
2257 :
2258 : LIBSAKURA_SYMBOL(Status) status;
2259 :
2260 20 : std::vector<float> weight(WeightIndex_kNum);
2261 10 : size_t const var_index = (sigma_value == "stddev") ? WeightIndex_kStddev : WeightIndex_kRms;
2262 :
2263 42 : for (vi->originChunks(); vi->moreChunks(); vi->nextChunk()) {
2264 64 : for (vi->origin(); vi->more(); vi->next()) {
2265 64 : Vector<Int> scans = vb->scan();
2266 64 : Vector<Double> times = vb->time();
2267 64 : Vector<Double> intervals = vb->timeInterval();
2268 64 : Vector<Int> beams = vb->feed1();
2269 64 : Vector<Int> antennas = vb->antenna1();
2270 64 : Vector<Int> data_spw = vb->spectralWindows();
2271 32 : size_t const num_chan = static_cast<size_t>(vb->nChannels());
2272 32 : size_t const num_pol = static_cast<size_t>(vb->nCorrelations());
2273 32 : size_t const num_row = static_cast<size_t>(vb->nRows());
2274 64 : Cube<Float> data_chunk(num_pol, num_chan, num_row);
2275 64 : SakuraAlignedArray<float> spec(num_chan);
2276 32 : float *spec_data = spec.data();
2277 64 : Cube<Bool> flag_chunk(num_pol, num_chan, num_row);
2278 64 : SakuraAlignedArray<bool> flag(num_chan);
2279 32 : bool *flag_data = flag.data();
2280 64 : SakuraAlignedArray<bool> mask(num_chan);
2281 32 : bool *mask_data = mask.data();
2282 64 : Matrix<Float> weight_matrix(num_pol, num_row, Array<Float>::uninitialized);
2283 :
2284 32 : bool new_nchan = false;
2285 32 : get_nchan_and_mask(recspw, data_spw, recchan, num_chan, nchan, in_mask, nchan_set, new_nchan);
2286 32 : if (new_nchan) {
2287 10 : map<size_t const, uint16_t>::iterator iter = max_orders.begin();
2288 30 : while (iter != max_orders.end()) {
2289 20 : get_baseline_context((*iter).first, (*iter).second,
2290 : num_chan, nchan, nchan_set,
2291 20 : ctx_indices, context_reservoir[(*iter).first]);
2292 20 : ++iter;
2293 : }
2294 : }
2295 : // get data/flag cubes (npol*nchan*nrow) from VisBuffer
2296 32 : get_data_cube_float(*vb, data_chunk);
2297 32 : get_flag_cube(*vb, flag_chunk);
2298 :
2299 : // get weight matrix (npol*nrow) from VisBuffer
2300 32 : if (update_weight) {
2301 8 : get_weight_matrix(*vb, weight_matrix);
2302 : }
2303 :
2304 : // loop over MS rows
2305 64 : for (size_t irow = 0; irow < num_row; ++irow) {
2306 32 : size_t idx = 0;
2307 72 : for (size_t ispw = 0; ispw < recspw.nelements(); ++ispw) {
2308 72 : if (data_spw[irow] == recspw[ispw]) {
2309 32 : idx = ispw;
2310 32 : break;
2311 : }
2312 : }
2313 :
2314 : // find a baseline table row (index) corresponding to this MS row
2315 : size_t idx_fit_param;
2316 32 : if (!parser.GetFitParameterIdx(times[irow], intervals[irow],
2317 32 : scans[irow], beams[irow], antennas[irow], data_spw[irow],
2318 : idx_fit_param)) {
2319 3 : for (size_t ipol = 0; ipol < num_pol; ++ipol) {
2320 2 : flag_spectrum_in_cube(flag_chunk, irow, ipol); //flag
2321 : }
2322 1 : continue;
2323 : }
2324 :
2325 : // loop over polarization
2326 93 : for (size_t ipol = 0; ipol < num_pol; ++ipol) {
2327 : bool apply;
2328 62 : Vector<double> coeff;
2329 62 : Vector<size_t> boundary;
2330 62 : std::vector<bool> mask_bltable;
2331 124 : BLParameterSet fit_param;
2332 62 : parser.GetFitParameterByIdx(idx_fit_param, ipol, apply, coeff, boundary,
2333 : mask_bltable, fit_param);
2334 62 : if (!apply) {
2335 1 : flag_spectrum_in_cube(flag_chunk, irow, ipol); //flag
2336 1 : continue;
2337 : }
2338 :
2339 : // get a channel mask from data cube
2340 : // (note that the variable 'mask' is flag in the next line
2341 : // actually, then it will be converted to real mask when
2342 : // taking AND with user-given mask info. this is just for
2343 : // saving memory usage...)
2344 61 : get_flag_from_cube(flag_chunk, irow, ipol, num_chan, flag_data);
2345 :
2346 : // skip spectrum if all channels flagged
2347 61 : if (allchannels_flagged(num_chan, flag_data)) {
2348 1 : continue;
2349 : }
2350 :
2351 : // get a spectrum from data cube
2352 60 : get_spectrum_from_cube(data_chunk, irow, ipol, num_chan, spec_data);
2353 :
2354 : // actual execution of single spectrum
2355 : map<size_t const, std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> >::iterator
2356 60 : iter = context_reservoir.find(fit_param.baseline_type);
2357 60 : if (iter == context_reservoir.end())
2358 0 : throw(AipsError("Invalid baseline type detected!"));
2359 : LIBSAKURA_SYMBOL(LSQFitContextFloat) * context =
2360 60 : (*iter).second[ctx_indices[idx]];
2361 : //cout << "Got context for type " << (*iter).first << ": idx=" << ctx_indices[idx] << endl;
2362 :
2363 120 : SakuraAlignedArray<double> coeff_storage(coeff);
2364 60 : double *coeff_data = coeff_storage.data();
2365 120 : SakuraAlignedArray<size_t> boundary_storage(boundary);
2366 60 : size_t *boundary_data = boundary_storage.data();
2367 120 : string subtract_funcname;
2368 60 : switch (static_cast<size_t>(fit_param.baseline_type)) {
2369 50 : case BaselineType_kPolynomial:
2370 : case BaselineType_kChebyshev:
2371 50 : status = LIBSAKURA_SYMBOL(SubtractPolynomialFloat)(
2372 : context, num_chan, spec_data, coeff.size(), coeff_data, spec_data);
2373 50 : subtract_funcname = "sakura_SubtractPolynomialFloat";
2374 50 : break;
2375 10 : case BaselineType_kCubicSpline:
2376 10 : status = LIBSAKURA_SYMBOL(SubtractCubicSplineFloat)(
2377 10 : context, num_chan, spec_data, boundary.size()-1,
2378 : reinterpret_cast<double (*)[4]>(coeff_data),
2379 : boundary_data, spec_data);
2380 10 : subtract_funcname = "sakura_SubtractCubicSplineFloat";
2381 10 : break;
2382 0 : default:
2383 0 : throw(AipsError("Unsupported baseline type."));
2384 : }
2385 60 : check_sakura_status(subtract_funcname, status);
2386 :
2387 : // set back a spectrum to data cube
2388 60 : set_spectrum_to_cube(data_chunk, irow, ipol, num_chan, spec_data);
2389 :
2390 60 : if (update_weight) {
2391 : // convert flag to mask by taking logical NOT of flag
2392 : // and then operate logical AND with in_mask and with mask from bltable
2393 131088 : for (size_t ichan = 0; ichan < num_chan; ++ichan) {
2394 131072 : mask_data[ichan] = in_mask[idx][ichan] && (!(flag_data[ichan])) && mask_bltable[ichan];
2395 : }
2396 16 : compute_weight(num_chan, spec_data, mask_data, weight);
2397 16 : set_weight_to_matrix(weight_matrix, irow, ipol, weight.at(var_index));
2398 : }
2399 : } // end of polarization loop
2400 :
2401 : } // end of chunk row loop
2402 : // write back data and flag cube to VisBuffer
2403 32 : if (update_weight) {
2404 8 : sdh_->fillCubeToOutputMs(vb, data_chunk, &flag_chunk, &weight_matrix);
2405 : } else {
2406 24 : sdh_->fillCubeToOutputMs(vb, data_chunk, &flag_chunk);
2407 : }
2408 : } // end of vi loop
2409 : } // end of chunk loop
2410 :
2411 10 : finalize_process();
2412 : // destroy baseline contexts
2413 10 : map<size_t const, std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> >::iterator ctxiter = context_reservoir.begin();
2414 30 : while (ctxiter != context_reservoir.end()) {
2415 20 : destroy_baseline_contexts (context_reservoir[(*ctxiter).first]);
2416 20 : ++ctxiter;
2417 : }
2418 10 : }
2419 :
2420 : // Fit line profile
2421 50 : void SingleDishMS::fitLine(string const& in_column_name,
2422 : string const& in_spw,
2423 : string const& /* in_pol */,
2424 : string const& fitfunc,
2425 : string const& in_nfit,
2426 : bool const linefinding,
2427 : float const threshold,
2428 : int const avg_limit,
2429 : int const minwidth,
2430 : vector<int> const& edge,
2431 : string const& tempfile_name,
2432 : string const& temp_out_ms_name) {
2433 :
2434 : // in_column = [FLOAT_DATA|DATA|CORRECTED_DATA]
2435 : // no iteration is necessary for the processing.
2436 : // procedure
2437 : // 1. iterate over MS
2438 : // 2. pick single spectrum from in_column (this is not actually necessary for simple scaling but for exibision purpose)
2439 : // 3. fit Gaussian or Lorentzian profile to each spectrum
2440 : // 4. write fitting results to outfile
2441 :
2442 150 : LogIO os(_ORIGIN);
2443 50 : os << "Fitting line profile with " << fitfunc << LogIO::POST;
2444 :
2445 50 : if (file_exists(temp_out_ms_name)) {
2446 0 : throw(AipsError("temporary ms file unexpectedly exists."));
2447 : }
2448 50 : if (file_exists(tempfile_name)) {
2449 0 : throw(AipsError("temporary file unexpectedly exists."));
2450 : }
2451 :
2452 100 : Block<Int> columns(1);
2453 50 : columns[0] = MS::DATA_DESC_ID;
2454 50 : prepare_for_process(in_column_name, temp_out_ms_name, columns, false);
2455 50 : vi::VisibilityIterator2 *vi = sdh_->getVisIter();
2456 50 : vi->setRowBlocking(kNRowBlocking);
2457 50 : vi::VisBuffer2 *vb = vi->getVisBuffer();
2458 100 : ofstream ofs(tempfile_name);
2459 :
2460 100 : Vector<Int> recspw;
2461 100 : Matrix<Int> recchan;
2462 100 : Vector<size_t> nchan;
2463 100 : Vector<Vector<Bool> > in_mask;
2464 100 : Vector<bool> nchan_set;
2465 50 : parse_spw(in_spw, recspw, recchan, nchan, in_mask, nchan_set);
2466 :
2467 100 : std::vector<size_t> nfit;
2468 50 : if (linefinding) {
2469 11 : os << "Defining line ranges using line finder. nfit will be ignored." << LogIO::POST;
2470 : } else {
2471 78 : std::vector<string> nfit_s = split_string(in_nfit, ',');
2472 39 : nfit.resize(nfit_s.size());
2473 82 : for (size_t i = 0; i < nfit_s.size(); ++i) {
2474 43 : nfit[i] = std::stoi(nfit_s[i]);
2475 : }
2476 : }
2477 :
2478 50 : size_t num_spec = 0;
2479 50 : size_t num_noline = 0;
2480 135 : for (vi->originChunks(); vi->moreChunks(); vi->nextChunk()) {
2481 170 : for (vi->origin(); vi->more(); vi->next()) {
2482 170 : Vector<Int> scans = vb->scan();
2483 170 : Vector<Double> times = vb->time();
2484 170 : Vector<Int> beams = vb->feed1();
2485 170 : Vector<Int> antennas = vb->antenna1();
2486 :
2487 170 : Vector<Int> data_spw = vb->spectralWindows();
2488 85 : size_t const num_chan = static_cast<size_t>(vb->nChannels());
2489 85 : size_t const num_pol = static_cast<size_t>(vb->nCorrelations());
2490 85 : size_t const num_row = static_cast<size_t>(vb->nRows());
2491 170 : Cube<Float> data_chunk(num_pol, num_chan, num_row);
2492 170 : SakuraAlignedArray<float> spec(num_chan);
2493 170 : Cube<Bool> flag_chunk(num_pol, num_chan, num_row);
2494 170 : SakuraAlignedArray<bool> mask(num_chan);
2495 : // CAUTION!!!
2496 : // data() method must be used with special care!!!
2497 85 : float *spec_data = spec.data();
2498 85 : bool *mask_data = mask.data();
2499 :
2500 85 : bool new_nchan = false;
2501 85 : get_nchan_and_mask(recspw, data_spw, recchan, num_chan, nchan, in_mask,
2502 : nchan_set, new_nchan);
2503 :
2504 : // get data/flag cubes (npol*nchan*nrow) from VisBuffer
2505 85 : get_data_cube_float(*vb, data_chunk);
2506 85 : get_flag_cube(*vb, flag_chunk);
2507 :
2508 : // loop over MS rows
2509 176 : for (size_t irow = 0; irow < num_row; ++irow) {
2510 91 : size_t idx = 0;
2511 441 : for (size_t ispw = 0; ispw < recspw.nelements(); ++ispw) {
2512 441 : if (data_spw[irow] == recspw[ispw]) {
2513 91 : idx = ispw;
2514 91 : break;
2515 : }
2516 : }
2517 :
2518 182 : std::vector<size_t> fitrange_start;
2519 91 : fitrange_start.clear();
2520 182 : std::vector<size_t> fitrange_end;
2521 91 : fitrange_end.clear();
2522 1100 : for (size_t i = 0; i < recchan.nrow(); ++i) {
2523 1009 : if (recchan.row(i)(0) == data_spw[irow]) {
2524 95 : fitrange_start.push_back(recchan.row(i)(1));
2525 95 : fitrange_end.push_back(recchan.row(i)(2));
2526 : }
2527 : }
2528 91 : if (!linefinding && nfit.size() != fitrange_start.size()) {
2529 : throw(AipsError(
2530 0 : "the number of elements of nfit and fitranges specified in spw must be identical."));
2531 : }
2532 :
2533 : // loop over polarization
2534 229 : for (size_t ipol = 0; ipol < num_pol; ++ipol) {
2535 : // get a channel mask from data cube
2536 : // (note that the variable 'mask' is flag in the next line
2537 : // actually, then it will be converted to real mask when
2538 : // taking AND with user-given mask info. this is just for
2539 : // saving memory usage...)
2540 138 : get_flag_from_cube(flag_chunk, irow, ipol, num_chan, mask_data);
2541 : // skip spectrum if all channels flagged
2542 138 : if (allchannels_flagged(num_chan, mask_data)) {
2543 0 : continue;
2544 : }
2545 138 : ++num_spec;
2546 :
2547 : // convert flag to mask by taking logical NOT of flag
2548 : // and then operate logical AND with in_mask
2549 533898 : for (size_t ichan = 0; ichan < num_chan; ++ichan) {
2550 533760 : mask_data[ichan] = in_mask[idx][ichan] && (!(mask_data[ichan]));
2551 : }
2552 : // get a spectrum from data cube
2553 138 : get_spectrum_from_cube(data_chunk, irow, ipol, num_chan, spec_data);
2554 :
2555 : // line finding. get fit mask (invert=false)
2556 138 : if (linefinding) {
2557 : list<pair<size_t, size_t>> line_ranges
2558 : = findLineAndGetRanges(num_chan, spec_data, mask_data,
2559 : threshold, avg_limit, minwidth,
2560 20 : edge, false);
2561 20 : if (line_ranges.size()==0) {
2562 0 : ++num_noline;
2563 0 : continue;
2564 : }
2565 20 : size_t nline = line_ranges.size();
2566 20 : fitrange_start.resize(nline);
2567 20 : fitrange_end.resize(nline);
2568 20 : nfit.resize(nline);
2569 20 : auto range=line_ranges.begin();
2570 50 : for (size_t iline=0; iline<nline; ++iline){
2571 30 : fitrange_start[iline] = (*range).first;
2572 30 : fitrange_end[iline] = (*range).second;
2573 30 : nfit[iline] = 1;
2574 30 : ++range;
2575 : }
2576 : }
2577 :
2578 276 : Vector<Float> x_;
2579 138 : x_.resize(num_chan);
2580 276 : Vector<Float> y_;
2581 138 : y_.resize(num_chan);
2582 276 : Vector<Bool> m_;
2583 138 : m_.resize(num_chan);
2584 533898 : for (size_t ichan = 0; ichan < num_chan; ++ichan) {
2585 533760 : x_[ichan] = static_cast<Float>(ichan);
2586 533760 : y_[ichan] = spec_data[ichan];
2587 : }
2588 276 : Vector<Float> parameters_;
2589 276 : Vector<Float> error_;
2590 :
2591 276 : PtrBlock<Function<Float>*> funcs_;
2592 276 : std::vector<std::string> funcnames_;
2593 276 : std::vector<int> funccomponents_;
2594 276 : std::string expr;
2595 138 : if (fitfunc == "gaussian") {
2596 118 : expr = "gauss";
2597 20 : } else if (fitfunc == "lorentzian") {
2598 20 : expr = "lorentz";
2599 : }
2600 :
2601 138 : bool any_converged = false;
2602 294 : for (size_t ifit = 0; ifit < nfit.size(); ++ifit) {
2603 156 : if (nfit[ifit] == 0)
2604 16 : continue;
2605 :
2606 140 : if (0 < ifit)
2607 18 : ofs << ":";
2608 :
2609 : //extract spec/mask within fitrange
2610 518028 : for (size_t ichan = 0; ichan < num_chan; ++ichan) {
2611 517888 : if ((fitrange_start[ifit] <= ichan)
2612 517888 : && (ichan <= fitrange_end[ifit])) {
2613 343353 : m_[ichan] = mask_data[ichan];
2614 : } else {
2615 174535 : m_[ichan] = false;
2616 : }
2617 : }
2618 :
2619 : //initial guesss
2620 280 : Vector<Float> peak;
2621 280 : Vector<Float> cent;
2622 280 : Vector<Float> fwhm;
2623 140 : peak.resize(nfit[ifit]);
2624 140 : cent.resize(nfit[ifit]);
2625 140 : fwhm.resize(nfit[ifit]);
2626 140 : if (nfit[ifit] == 1) {
2627 132 : Float sum = 0.0;
2628 132 : Float max_spec = fabs(y_[fitrange_start[ifit]]);
2629 132 : Float max_spec_x = x_[fitrange_start[ifit]];
2630 132 : bool is_positive = true;
2631 277949 : for (size_t ichan = fitrange_start[ifit];
2632 277949 : ichan <= fitrange_end[ifit]; ++ichan) {
2633 277817 : sum += y_[ichan];
2634 277817 : if (max_spec < fabs(y_[ichan])) {
2635 6769 : max_spec = fabs(y_[ichan]);
2636 6769 : max_spec_x = x_[ichan];
2637 6769 : is_positive = (fabs(y_[ichan]) == y_[ichan]);
2638 : }
2639 : }
2640 132 : peak[0] = max_spec * (is_positive ? 1 : -1);
2641 132 : cent[0] = max_spec_x;
2642 132 : fwhm[0] = fabs(sum / max_spec * 0.7);
2643 : } else {
2644 8 : size_t x_start = fitrange_start[ifit];
2645 8 : size_t x_width = (fitrange_end[ifit] - fitrange_start[ifit]) / nfit[ifit];
2646 8 : size_t x_end = x_start + x_width;
2647 24 : for (size_t icomp = 0; icomp < nfit[ifit]; ++icomp) {
2648 16 : if (icomp == nfit[ifit] - 1) {
2649 8 : x_end = fitrange_end[ifit] + 1;
2650 : }
2651 :
2652 16 : Float sum = 0.0;
2653 16 : Float max_spec = fabs(y_[x_start]);
2654 16 : Float max_spec_x = x_[x_start];
2655 16 : bool is_positive = true;
2656 65552 : for (size_t ichan = x_start; ichan < x_end; ++ichan) {
2657 65536 : sum += y_[ichan];
2658 65536 : if (max_spec < fabs(y_[ichan])) {
2659 878 : max_spec = fabs(y_[ichan]);
2660 878 : max_spec_x = x_[ichan];
2661 878 : is_positive = (fabs(y_[ichan]) == y_[ichan]);
2662 : }
2663 : }
2664 16 : peak[icomp] = max_spec * (is_positive ? 1 : -1);
2665 16 : cent[icomp] = max_spec_x;
2666 16 : fwhm[icomp] = fabs(sum / max_spec * 0.7);
2667 :
2668 16 : x_start += x_width;
2669 16 : x_end += x_width;
2670 : }
2671 : }
2672 :
2673 : //fitter setup
2674 140 : funcs_.resize(nfit[ifit]);
2675 140 : funcnames_.clear();
2676 140 : funccomponents_.clear();
2677 288 : for (size_t icomp = 0; icomp < funcs_.nelements(); ++icomp) {
2678 148 : if (expr == "gauss") {
2679 120 : funcs_[icomp] = new Gaussian1D<Float>();
2680 28 : } else if (expr == "lorentz") {
2681 28 : funcs_[icomp] = new Lorentzian1D<Float>();
2682 : }
2683 148 : (funcs_[icomp]->parameters())[0] = peak[icomp]; //initial guess (peak)
2684 148 : (funcs_[icomp]->parameters())[1] = cent[icomp]; //initial guess (centre)
2685 148 : (funcs_[icomp]->parameters())[2] = fwhm[icomp]; //initial guess (fwhm)
2686 148 : funcnames_.push_back(expr);
2687 148 : funccomponents_.push_back(3);
2688 : }
2689 :
2690 : //actual fitting
2691 280 : NonLinearFitLM<Float> fitter;
2692 280 : CompoundFunction<Float> func;
2693 288 : for (size_t icomp = 0; icomp < funcs_.nelements(); ++icomp) {
2694 148 : func.addFunction(*funcs_[icomp]);
2695 : }
2696 140 : fitter.setFunction(func);
2697 140 : fitter.setMaxIter(50 + 10 * funcs_.nelements());
2698 140 : fitter.setCriteria(0.001); // Convergence criterium
2699 :
2700 140 : parameters_.resize();
2701 140 : parameters_ = fitter.fit(x_, y_, &m_);
2702 140 : any_converged |= fitter.converged();
2703 : // if (!fitter.converged()) {
2704 : // throw(AipsError("Failed in fitting. Fitter did not converge."));
2705 : // }
2706 140 : error_.resize();
2707 140 : error_ = fitter.errors();
2708 :
2709 : //write best-fit parameters to tempfile/outfile
2710 288 : for (size_t icomp = 0; icomp < funcs_.nelements(); ++icomp) {
2711 148 : if (0 < icomp)
2712 8 : ofs << ":";
2713 148 : size_t offset = 3 * icomp;
2714 148 : ofs.precision(4);
2715 148 : ofs.setf(ios::fixed);
2716 148 : ofs << scans[irow] << "," // scanID
2717 148 : << times[irow] << "," // time
2718 148 : << antennas[irow] << "," // antennaID
2719 148 : << beams[irow] << "," // beamID
2720 148 : << data_spw[irow] << "," // spwID
2721 148 : << ipol << ","; // polID
2722 148 : ofs.precision(8);
2723 148 : ofs << parameters_[offset + 1] << "," << error_[offset + 1] << "," // cent
2724 148 : << parameters_[offset + 0] << "," << error_[offset + 0] << "," // peak
2725 148 : << parameters_[offset + 2] << "," << error_[offset + 2]; // fwhm
2726 : }
2727 : } //end of nfit loop
2728 138 : ofs << "\n";
2729 : // count up spectra w/o any line fit
2730 138 : if (!any_converged) ++num_noline;
2731 :
2732 : } //end of polarization loop
2733 : } // end of MS row loop
2734 : } //end of vi loop
2735 : } //end of chunk loop
2736 :
2737 50 : if (num_noline==num_spec) {
2738 : os << LogIO::WARN
2739 2 : << "Fitter did not converge on any fit components." << LogIO::POST;
2740 : }
2741 48 : else if (num_noline > 0) {
2742 : os << "No convergence for fitting to " << num_noline
2743 0 : << " out of " << num_spec << " spectra" << LogIO::POST;
2744 : }
2745 :
2746 50 : finalize_process();
2747 50 : ofs.close();
2748 50 : }
2749 :
2750 : //Subtract baseline by per spectrum fitting parameters
2751 81 : void SingleDishMS::subtractBaselineVariable(string const& in_column_name,
2752 : string const& out_ms_name,
2753 : string const& out_bloutput_name,
2754 : bool const& do_subtract,
2755 : string const& in_spw,
2756 : bool const& update_weight,
2757 : string const& sigma_value,
2758 : string const& param_file,
2759 : bool const& verbose) {
2760 :
2761 243 : LogIO os(_ORIGIN);
2762 : os << "Fitting and subtracting baseline using parameters in file "
2763 81 : << param_file << LogIO::POST;
2764 :
2765 162 : Block<Int> columns(1);
2766 81 : columns[0] = MS::DATA_DESC_ID;
2767 81 : prepare_for_process(in_column_name, out_ms_name, columns, false);
2768 81 : vi::VisibilityIterator2 *vi = sdh_->getVisIter();
2769 81 : vi->setRowBlocking(kNRowBlocking);
2770 81 : vi::VisBuffer2 *vb = vi->getVisBuffer();
2771 :
2772 81 : split_bloutputname(out_bloutput_name);
2773 81 : bool write_baseline_csv = (bloutputname_csv != "");
2774 81 : bool write_baseline_text = (bloutputname_text != "");
2775 81 : bool write_baseline_table = (bloutputname_table != "");
2776 162 : ofstream ofs_csv;
2777 162 : ofstream ofs_txt;
2778 81 : BaselineTable *bt = 0;
2779 :
2780 81 : if (write_baseline_csv) {
2781 27 : ofs_csv.open(bloutputname_csv.c_str());
2782 : }
2783 81 : if (write_baseline_text) {
2784 12 : ofs_txt.open(bloutputname_text.c_str(), std::ios::app);
2785 : }
2786 81 : if (write_baseline_table) {
2787 21 : bt = new BaselineTable(vi->ms());
2788 : }
2789 :
2790 162 : Vector<Int> recspw;
2791 162 : Matrix<Int> recchan;
2792 162 : Vector<size_t> nchan;
2793 162 : Vector<Vector<Bool> > in_mask;
2794 162 : Vector<bool> nchan_set;
2795 81 : parse_spw(in_spw, recspw, recchan, nchan, in_mask, nchan_set);
2796 :
2797 : // parse fitting parameters in the text file
2798 162 : BLParameterParser parser(param_file);
2799 162 : std::vector<size_t> baseline_types = parser.get_function_types();
2800 162 : map<size_t const, uint16_t> max_orders;
2801 289 : for (size_t i = 0; i < baseline_types.size(); ++i) {
2802 208 : max_orders[baseline_types[i]] = parser.get_max_order(baseline_types[i]);
2803 : }
2804 : { //DEBUG ouput
2805 81 : os << LogIO::DEBUG1 << "Baseline Types = " << baseline_types << LogIO::POST;
2806 81 : os << LogIO::DEBUG1 << "Max Orders:" << LogIO::POST;
2807 81 : map<size_t const, uint16_t>::iterator iter = max_orders.begin();
2808 289 : while (iter != max_orders.end()) {
2809 416 : os << LogIO::DEBUG1 << "- type " << (*iter).first << ": "
2810 208 : << (*iter).second << LogIO::POST;
2811 208 : ++iter;
2812 : }
2813 : }
2814 :
2815 : // Baseline Contexts reservoir
2816 : // key: Sakura_BaselineType enum,
2817 : // value: a vector of Sakura_BaselineContextFloat for various nchans
2818 162 : map<size_t const, std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> > context_reservoir;
2819 : {
2820 81 : map<size_t const, uint16_t>::iterator iter = max_orders.begin();
2821 289 : while (iter != max_orders.end()) {
2822 208 : context_reservoir[(*iter).first] = std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *>();
2823 208 : ++iter;
2824 : }
2825 : }
2826 :
2827 162 : Vector<size_t> ctx_indices(recspw.size(), 0ul);
2828 : //stores the number of channels of corresponding elements in contexts list.
2829 : // WORKAROUND for absense of the way to get num_bases_data in context.
2830 162 : vector<size_t> ctx_nchans;
2831 :
2832 : LIBSAKURA_SYMBOL(Status) status;
2833 : LIBSAKURA_SYMBOL(LSQFitStatus) bl_status;
2834 :
2835 162 : std::vector<float> weight(WeightIndex_kNum);
2836 81 : size_t const var_index = (sigma_value == "stddev") ? WeightIndex_kStddev : WeightIndex_kRms;
2837 :
2838 300 : for (vi->originChunks(); vi->moreChunks(); vi->nextChunk()) {
2839 438 : for (vi->origin(); vi->more(); vi->next()) {
2840 438 : Vector<Int> scans = vb->scan();
2841 438 : Vector<Double> times = vb->time();
2842 438 : Vector<Int> beams = vb->feed1();
2843 438 : Vector<Int> antennas = vb->antenna1();
2844 438 : Vector<Int> data_spw = vb->spectralWindows();
2845 219 : size_t const num_chan = static_cast<size_t>(vb->nChannels());
2846 219 : size_t const num_pol = static_cast<size_t>(vb->nCorrelations());
2847 219 : size_t const num_row = static_cast<size_t>(vb->nRows());
2848 438 : auto orig_rows = vb->rowIds();
2849 438 : Cube<Float> data_chunk(num_pol, num_chan, num_row);
2850 438 : SakuraAlignedArray<float> spec(num_chan);
2851 438 : Cube<Bool> flag_chunk(num_pol, num_chan, num_row);
2852 438 : SakuraAlignedArray<bool> flag(num_chan);
2853 438 : SakuraAlignedArray<bool> mask(num_chan);
2854 438 : SakuraAlignedArray<bool> mask_after_clipping(num_chan);
2855 : // CAUTION!!!
2856 : // data() method must be used with special care!!!
2857 219 : float *spec_data = spec.data();
2858 219 : bool *flag_data = flag.data();
2859 219 : bool *mask_data = mask.data();
2860 219 : bool *mask_after_clipping_data = mask_after_clipping.data();
2861 438 : Matrix<Float> weight_matrix(num_pol, num_row, Array<Float>::uninitialized);
2862 :
2863 219 : uInt final_mask[num_pol];
2864 219 : uInt final_mask_after_clipping[num_pol];
2865 219 : final_mask[0] = 0;
2866 219 : final_mask[1] = 0;
2867 219 : final_mask_after_clipping[0] = 0;
2868 219 : final_mask_after_clipping[1] = 0;
2869 :
2870 219 : bool new_nchan = false;
2871 219 : get_nchan_and_mask(recspw, data_spw, recchan, num_chan, nchan, in_mask, nchan_set, new_nchan);
2872 : // check if context should be created once per chunk
2873 : // in the first actual excution of baseline.
2874 219 : bool check_context = true;
2875 :
2876 : // get data/flag cubes (npol*nchan*nrow) from VisBuffer
2877 219 : get_data_cube_float(*vb, data_chunk);
2878 219 : get_flag_cube(*vb, flag_chunk);
2879 :
2880 : // get weight matrix (npol*nrow) from VisBuffer
2881 219 : if (update_weight) {
2882 6 : get_weight_matrix(*vb, weight_matrix);
2883 : }
2884 :
2885 : // loop over MS rows
2886 2814 : for (size_t irow = 0; irow < num_row; ++irow) {
2887 2595 : size_t idx = 0;
2888 2856 : for (size_t ispw = 0; ispw < recspw.nelements(); ++ispw) {
2889 2856 : if (data_spw[irow] == recspw[ispw]) {
2890 2595 : idx = ispw;
2891 2595 : break;
2892 : }
2893 : }
2894 :
2895 : //prepare varables for writing baseline table
2896 2595 : Array<Bool> apply_mtx(IPosition(2, num_pol, 1), true);
2897 2595 : Array<uInt> bltype_mtx(IPosition(2, num_pol, 1));
2898 2595 : Array<Int> fpar_mtx(IPosition(2, num_pol, 1));
2899 2595 : std::vector<std::vector<double> > ffpar_mtx_tmp(num_pol);
2900 2595 : std::vector<std::vector<uInt> > masklist_mtx_tmp(num_pol);
2901 2595 : std::vector<std::vector<double> > coeff_mtx_tmp(num_pol);
2902 2595 : Array<Float> rms_mtx(IPosition(2, num_pol, 1));
2903 2595 : Array<Float> cthres_mtx(IPosition(2, num_pol, 1));
2904 2595 : Array<uInt> citer_mtx(IPosition(2, num_pol, 1));
2905 2595 : Array<Bool> uself_mtx(IPosition(2, num_pol, 1));
2906 2595 : Array<Float> lfthres_mtx(IPosition(2, num_pol, 1));
2907 2595 : Array<uInt> lfavg_mtx(IPosition(2, num_pol, 1));
2908 2595 : Array<uInt> lfedge_mtx(IPosition(2, num_pol, 2));
2909 :
2910 2595 : size_t num_apply_true = 0;
2911 2595 : size_t num_ffpar_max = 0;
2912 2595 : size_t num_masklist_max = 0;
2913 2595 : size_t num_coeff_max = 0;
2914 2595 : std::vector<size_t> numcoeff(num_pol);
2915 :
2916 : // loop over polarization
2917 5377 : for (size_t ipol = 0; ipol < num_pol; ++ipol) {
2918 : // get a channel mask from data cube
2919 : // (note that the variable 'mask' is flag in the next line
2920 : // actually, then it will be converted to real mask when
2921 : // taking AND with user-given mask info. this is just for
2922 : // saving memory usage...)
2923 2782 : get_flag_from_cube(flag_chunk, irow, ipol, num_chan, flag_data);
2924 : // skip spectrum if all channels flagged
2925 2782 : if (allchannels_flagged(num_chan, flag_data)) {
2926 1736 : os << LogIO::DEBUG1 << "Row " << orig_rows[irow] << ", Pol " << ipol
2927 868 : << ": All channels flagged. Skipping." << LogIO::POST;
2928 868 : apply_mtx[0][ipol] = false;
2929 2541 : continue;
2930 : }
2931 :
2932 : // convert flag to mask by taking logical NOT of flag
2933 : // and then operate logical AND with in_mask
2934 9365370 : for (size_t ichan = 0; ichan < num_chan; ++ichan) {
2935 9363456 : mask_data[ichan] = in_mask[idx][ichan] && (!(flag_data[ichan]));
2936 : }
2937 : // get fitting parameter
2938 3828 : BLParameterSet fit_param;
2939 1914 : if (!parser.GetFitParameter(orig_rows[irow], ipol, fit_param)) { //no fit requrested
2940 1673 : flag_spectrum_in_cube(flag_chunk, irow, ipol);
2941 3346 : os << LogIO::DEBUG1 << "Row " << orig_rows[irow] << ", Pol " << ipol
2942 1673 : << ": Fit not requested. Skipping." << LogIO::POST;
2943 1673 : apply_mtx[0][ipol] = false;
2944 1673 : continue;
2945 : }
2946 241 : if (verbose) {
2947 0 : os << "Fitting Parameter" << LogIO::POST;
2948 0 : os << "[ROW " << orig_rows[irow] << " (nchan " << num_chan << ")" << ", POL" << ipol << "]"
2949 0 : << LogIO::POST;
2950 0 : fit_param.PrintSummary();
2951 : }
2952 : // Create contexts when actually subtract baseine for the first time (if not yet exist)
2953 241 : if (check_context) {
2954 : // Generate context for all necessary baseline types
2955 128 : map<size_t const, uint16_t>::iterator iter = max_orders.begin();
2956 499 : while (iter != max_orders.end()) {
2957 371 : get_baseline_context((*iter).first, (*iter).second, num_chan, idx,
2958 371 : ctx_indices, ctx_nchans, context_reservoir[(*iter).first]);
2959 371 : ++iter;
2960 : }
2961 128 : check_context = false;
2962 : }
2963 : // get mask from BLParameterset and create composit mask
2964 241 : if (fit_param.baseline_mask != "") {
2965 154 : stringstream local_spw;
2966 77 : local_spw << data_spw[irow] << ":" << fit_param.baseline_mask;
2967 154 : Record selrec = sdh_->getSelRec(local_spw.str());
2968 154 : Matrix<Int> local_rec_chan = selrec.asArrayInt("channel");
2969 154 : Vector<Bool> local_mask(num_chan, false);
2970 77 : get_mask_from_rec(data_spw[irow], local_rec_chan, local_mask, false);
2971 630861 : for (size_t ichan = 0; ichan < num_chan; ++ichan) {
2972 630784 : mask_data[ichan] = mask_data[ichan] && local_mask[ichan];
2973 : }
2974 : }
2975 : // check for composit mask and flag if no valid channel to fit
2976 241 : if (NValidMask(num_chan, mask_data) == 0) {
2977 0 : flag_spectrum_in_cube(flag_chunk, irow, ipol);
2978 0 : apply_mtx[0][ipol] = false;
2979 0 : os << LogIO::DEBUG1 << "Row " << orig_rows[irow] << ", Pol " << ipol
2980 0 : << ": No valid channel to fit. Skipping" << LogIO::POST;
2981 0 : continue;
2982 : }
2983 : // get a spectrum from data cube
2984 241 : get_spectrum_from_cube(data_chunk, irow, ipol, num_chan, spec_data);
2985 :
2986 : // get baseline context
2987 : map<size_t const, std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> >::iterator
2988 241 : iter = context_reservoir.find(fit_param.baseline_type);
2989 241 : if (iter == context_reservoir.end()) {
2990 0 : throw(AipsError("Invalid baseline type detected!"));
2991 : }
2992 241 : LIBSAKURA_SYMBOL(LSQFitContextFloat) * context = (*iter).second[ctx_indices[idx]];
2993 :
2994 : // Number of coefficients to fit this spectrum
2995 : size_t num_coeff;
2996 241 : size_t bltype = static_cast<size_t>(fit_param.baseline_type);
2997 241 : switch (bltype) {
2998 177 : case BaselineType_kPolynomial:
2999 : case BaselineType_kChebyshev:
3000 177 : status = LIBSAKURA_SYMBOL(GetNumberOfCoefficientsFloat)(context, fit_param.order, &num_coeff);
3001 177 : check_sakura_status("sakura_GetNumberOfCoefficientsFloat", status);
3002 177 : break;
3003 64 : case BaselineType_kCubicSpline:
3004 64 : num_coeff = 4 * fit_param.npiece;
3005 64 : break;
3006 0 : default:
3007 0 : throw(AipsError("Unsupported baseline type."));
3008 : }
3009 241 : numcoeff[ipol] = num_coeff;
3010 :
3011 : // Final check of the valid number of channels
3012 241 : size_t num_min = (bltype == BaselineType_kCubicSpline) ? fit_param.npiece + 3 : num_coeff;
3013 241 : if (NValidMask(num_chan, mask_data) < num_min) {
3014 0 : flag_spectrum_in_cube(flag_chunk, irow, ipol);
3015 0 : apply_mtx[0][ipol] = false;
3016 : os << LogIO::WARN
3017 : << "Too few valid channels to fit. Skipping Antenna "
3018 0 : << antennas[irow] << ", Beam " << beams[irow] << ", SPW "
3019 0 : << data_spw[irow] << ", Pol " << ipol << ", Time "
3020 0 : << MVTime(times[irow] / 24. / 3600.).string(MVTime::YMD, 8)
3021 0 : << LogIO::POST;
3022 0 : continue;
3023 : }
3024 :
3025 : // actual execution of single spectrum
3026 : float rms;
3027 241 : size_t num_boundary = 0;
3028 241 : if (bltype == BaselineType_kCubicSpline) {
3029 64 : num_boundary = fit_param.npiece+1;
3030 : }
3031 482 : SakuraAlignedArray<size_t> boundary(num_boundary);
3032 241 : size_t *boundary_data = boundary.data();
3033 :
3034 241 : if (write_baseline_text || write_baseline_csv || write_baseline_table) {
3035 159 : num_apply_true++;
3036 159 : bltype_mtx[0][ipol] = (uInt)fit_param.baseline_type;
3037 : Int fpar_tmp;
3038 159 : switch (bltype) {
3039 115 : case BaselineType_kPolynomial:
3040 : case BaselineType_kChebyshev:
3041 115 : fpar_tmp = (Int)fit_param.order;
3042 115 : break;
3043 44 : case BaselineType_kCubicSpline:
3044 44 : fpar_tmp = (Int)fit_param.npiece;
3045 44 : break;
3046 0 : default:
3047 0 : throw(AipsError("Unsupported baseline type."));
3048 : }
3049 159 : fpar_mtx[0][ipol] = fpar_tmp;
3050 :
3051 159 : if (num_coeff > num_coeff_max) {
3052 89 : num_coeff_max = num_coeff;
3053 : }
3054 318 : SakuraAlignedArray<double> coeff(num_coeff);
3055 : // CAUTION!!!
3056 : // data() method must be used with special care!!!
3057 159 : double *coeff_data = coeff.data();
3058 318 : string get_coeff_funcname;
3059 159 : switch (bltype) {
3060 115 : case BaselineType_kPolynomial:
3061 : case BaselineType_kChebyshev:
3062 345 : status = LIBSAKURA_SYMBOL(LSQFitPolynomialFloat)(
3063 115 : context, fit_param.order,
3064 : num_chan, spec_data, mask_data,
3065 115 : fit_param.clip_threshold_sigma, fit_param.num_fitting_max,
3066 : num_coeff, coeff_data, nullptr, nullptr,
3067 : mask_after_clipping_data, &rms, &bl_status);
3068 :
3069 942195 : for (size_t i = 0; i < num_chan; ++i) {
3070 942080 : if (mask_data[i] == false) {
3071 210640 : final_mask[ipol] += 1;
3072 : }
3073 942080 : if (mask_after_clipping_data[i] == false) {
3074 216221 : final_mask_after_clipping[ipol] += 1;
3075 : }
3076 : }
3077 :
3078 115 : get_coeff_funcname = "sakura_LSQFitPolynomialFloat";
3079 115 : break;
3080 44 : case BaselineType_kCubicSpline:
3081 44 : status = LIBSAKURA_SYMBOL(LSQFitCubicSplineFloat)(
3082 : context, fit_param.npiece,
3083 : num_chan, spec_data, mask_data,
3084 44 : fit_param.clip_threshold_sigma, fit_param.num_fitting_max,
3085 : reinterpret_cast<double (*)[4]>(coeff_data), nullptr, nullptr,
3086 : mask_after_clipping_data, &rms, boundary_data, &bl_status);
3087 :
3088 360492 : for (size_t i = 0; i < num_chan; ++i) {
3089 360448 : if (mask_data[i] == false) {
3090 83990 : final_mask[ipol] += 1;
3091 : }
3092 360448 : if (mask_after_clipping_data[i] == false) {
3093 86619 : final_mask_after_clipping[ipol] += 1;
3094 : }
3095 : }
3096 :
3097 44 : get_coeff_funcname = "sakura_LSQFitCubicSplineFloat";
3098 44 : break;
3099 0 : default:
3100 0 : throw(AipsError("Unsupported baseline type."));
3101 : }
3102 159 : check_sakura_status(get_coeff_funcname, status);
3103 :
3104 159 : size_t num_ffpar = get_num_coeff_bloutput(fit_param.baseline_type, fit_param.npiece, num_ffpar_max);
3105 159 : ffpar_mtx_tmp[ipol].clear();
3106 250 : for (size_t ipiece = 0; ipiece < num_ffpar; ++ipiece) {
3107 91 : ffpar_mtx_tmp[ipol].push_back(boundary_data[ipiece]);
3108 : }
3109 :
3110 159 : coeff_mtx_tmp[ipol].clear();
3111 586 : for (size_t icoeff = 0; icoeff < num_coeff; ++icoeff) {
3112 427 : coeff_mtx_tmp[ipol].push_back(coeff_data[icoeff]);
3113 : }
3114 :
3115 318 : Vector<uInt> masklist;
3116 159 : get_masklist_from_mask(num_chan, mask_after_clipping_data, masklist);
3117 159 : if (masklist.size() > num_masklist_max) {
3118 87 : num_masklist_max = masklist.size();
3119 : }
3120 159 : masklist_mtx_tmp[ipol].clear();
3121 851 : for (size_t imask = 0; imask < masklist.size(); ++imask) {
3122 692 : masklist_mtx_tmp[ipol].push_back(masklist[imask]);
3123 : }
3124 :
3125 318 : string subtract_funcname;
3126 159 : switch (fit_param.baseline_type) {
3127 115 : case BaselineType_kPolynomial:
3128 : case BaselineType_kChebyshev:
3129 115 : status = LIBSAKURA_SYMBOL(SubtractPolynomialFloat)(
3130 : context, num_chan, spec_data, num_coeff, coeff_data,
3131 : spec_data);
3132 115 : subtract_funcname = "sakura_SubtractPolynomialFloat";
3133 115 : break;
3134 44 : case BaselineType_kCubicSpline:
3135 44 : status = LIBSAKURA_SYMBOL(SubtractCubicSplineFloat)(
3136 : context,
3137 : num_chan, spec_data, fit_param.npiece, reinterpret_cast<double (*)[4]>(coeff_data),
3138 : boundary_data, spec_data);
3139 44 : subtract_funcname = "sakura_SubtractCubicSplineFloat";
3140 44 : break;
3141 0 : default:
3142 0 : throw(AipsError("Unsupported baseline type."));
3143 : }
3144 159 : check_sakura_status(subtract_funcname, status);
3145 :
3146 159 : rms_mtx[0][ipol] = rms;
3147 :
3148 159 : cthres_mtx[0][ipol] = fit_param.clip_threshold_sigma;
3149 159 : citer_mtx[0][ipol] = (uInt)fit_param.num_fitting_max - 1;
3150 159 : uself_mtx[0][ipol] = (Bool)fit_param.line_finder.use_line_finder;
3151 159 : lfthres_mtx[0][ipol] = fit_param.line_finder.threshold;
3152 159 : lfavg_mtx[0][ipol] = fit_param.line_finder.chan_avg_limit;
3153 477 : for (size_t iedge = 0; iedge < 2; ++iedge) {
3154 318 : lfedge_mtx[iedge][ipol] = fit_param.line_finder.edge[iedge];
3155 159 : }
3156 :
3157 : } else {
3158 164 : string subtract_funcname;
3159 82 : switch (fit_param.baseline_type) {
3160 62 : case BaselineType_kPolynomial:
3161 : case BaselineType_kChebyshev:
3162 186 : status = LIBSAKURA_SYMBOL(LSQFitPolynomialFloat)(
3163 62 : context, fit_param.order,
3164 : num_chan, spec_data, mask_data,
3165 62 : fit_param.clip_threshold_sigma, fit_param.num_fitting_max,
3166 : num_coeff, nullptr, nullptr, spec_data,
3167 : mask_after_clipping_data, &rms, &bl_status);
3168 62 : subtract_funcname = "sakura_LSQFitPolynomialFloat";
3169 62 : break;
3170 20 : case BaselineType_kCubicSpline:
3171 20 : status = LIBSAKURA_SYMBOL(LSQFitCubicSplineFloat)(
3172 : context, fit_param.npiece,
3173 : num_chan, spec_data, mask_data,
3174 20 : fit_param.clip_threshold_sigma, fit_param.num_fitting_max,
3175 : nullptr, nullptr, spec_data,
3176 : mask_after_clipping_data, &rms, boundary_data, &bl_status);
3177 20 : subtract_funcname = "sakura_LSQFitCubicSplineFloat";
3178 20 : break;
3179 0 : default:
3180 0 : throw(AipsError("Unsupported baseline type."));
3181 : }
3182 82 : check_sakura_status(subtract_funcname, status);
3183 : }
3184 : // set back a spectrum to data cube
3185 241 : if (do_subtract) {
3186 233 : set_spectrum_to_cube(data_chunk, irow, ipol, num_chan, spec_data);
3187 : }
3188 :
3189 241 : if (update_weight) {
3190 12 : compute_weight(num_chan, spec_data, mask_data, weight);
3191 12 : set_weight_to_matrix(weight_matrix, irow, ipol, weight.at(var_index));
3192 : }
3193 : } // end of polarization loop
3194 :
3195 : // output results of fitting
3196 2595 : if (num_apply_true == 0) continue;
3197 :
3198 172 : Array<Float> ffpar_mtx(IPosition(2, num_pol, num_ffpar_max));
3199 86 : set_matrix_for_bltable<double, Float>(num_pol, num_ffpar_max,
3200 : ffpar_mtx_tmp, ffpar_mtx);
3201 172 : Array<uInt> masklist_mtx(IPosition(2, num_pol, num_masklist_max));
3202 86 : set_matrix_for_bltable<uInt, uInt>(num_pol, num_masklist_max,
3203 : masklist_mtx_tmp, masklist_mtx);
3204 172 : Array<Float> coeff_mtx(IPosition(2, num_pol, num_coeff_max));
3205 86 : set_matrix_for_bltable<double, Float>(num_pol, num_coeff_max,
3206 : coeff_mtx_tmp, coeff_mtx);
3207 172 : Matrix<uInt> masklist_mtx2 = masklist_mtx;
3208 172 : Matrix<Bool> apply_mtx2 = apply_mtx;
3209 :
3210 86 : if (write_baseline_table) {
3211 122 : bt->appenddata((uInt)scans[irow], (uInt)beams[irow], (uInt)antennas[irow],
3212 61 : (uInt)data_spw[irow], 0, times[irow],
3213 : apply_mtx, bltype_mtx, fpar_mtx, ffpar_mtx, masklist_mtx,
3214 : coeff_mtx, rms_mtx, (uInt)num_chan, cthres_mtx, citer_mtx,
3215 : uself_mtx, lfthres_mtx, lfavg_mtx, lfedge_mtx);
3216 : }
3217 :
3218 86 : if (write_baseline_text) {
3219 69 : for (size_t ipol = 0; ipol < num_pol; ++ipol) {
3220 46 : if (apply_mtx2(ipol, 0) == false) continue;
3221 :
3222 43 : ofs_txt << "Scan" << '[' << (uInt)scans[irow] << ']' << ' '
3223 43 : << "Beam" << '[' << (uInt)beams[irow] << ']' << ' '
3224 43 : << "Spw" << '[' << (uInt)data_spw[irow] << ']' << ' '
3225 43 : << "Pol" << '[' << ipol << ']' << ' '
3226 86 : << "Time" << '[' << MVTime(times[irow]/ 24. / 3600.).string(MVTime::YMD, 8) << ']'
3227 43 : << endl;
3228 43 : ofs_txt << endl;
3229 43 : ofs_txt << "Fitter range = " << '[';
3230 :
3231 248 : for (size_t imasklist = 0; imasklist < num_masklist_max / 2; ++imasklist) {
3232 205 : if (imasklist == 0) {
3233 43 : ofs_txt << '[' << masklist_mtx2(ipol, 2 * imasklist) << ';'
3234 43 : << masklist_mtx2(ipol, 2 * imasklist + 1) << ']';
3235 : }
3236 205 : if (imasklist >= 1
3237 317 : && (0 != masklist_mtx2(ipol, 2 * imasklist)
3238 112 : && 0 != masklist_mtx2(ipol, 2 * imasklist + 1))) {
3239 112 : ofs_txt << ",[" << masklist_mtx2(ipol, 2 * imasklist) << ','
3240 112 : << masklist_mtx2(ipol, 2 * imasklist + 1) << ']';
3241 : }
3242 : }
3243 :
3244 43 : ofs_txt << ']' << endl;
3245 43 : ofs_txt << endl;
3246 :
3247 129 : Matrix<uInt> bltype_mtx2 = bltype_mtx[0][ipol];
3248 129 : Matrix<Int> fpar_mtx2 = fpar_mtx[0][ipol];
3249 129 : Matrix<Float> rms_mtx2 = rms_mtx[0][ipol];
3250 86 : string bltype_name;
3251 86 : string blparam_name = "order";
3252 43 : if (bltype_mtx2(0, 0) == (uInt)0) {
3253 12 : bltype_name = "poly";
3254 31 : } else if (bltype_mtx2(0, 0) == (uInt)1) {
3255 17 : bltype_name = "chebyshev";
3256 14 : } else if (bltype_mtx2(0, 0) == (uInt)2) {
3257 14 : bltype_name = "cspline";
3258 14 : blparam_name = "npiece";
3259 : }
3260 :
3261 : ofs_txt << "Baseline parameters Function = "
3262 43 : << bltype_name << " " << blparam_name << " = "
3263 43 : << fpar_mtx2(0, 0) << endl;
3264 43 : ofs_txt << endl;
3265 43 : ofs_txt << "Results of baseline fit" << endl;
3266 43 : ofs_txt << endl;
3267 86 : Matrix<Float> coeff_mtx2 = coeff_mtx;
3268 186 : for (size_t icoeff = 0; icoeff < numcoeff[ipol]; ++icoeff) {
3269 143 : ofs_txt << "p" << icoeff << " = "
3270 143 : << setprecision(8) << coeff_mtx2(ipol, icoeff) << " ";
3271 : }
3272 43 : ofs_txt << endl;
3273 43 : ofs_txt << endl;
3274 43 : ofs_txt << "rms = ";
3275 43 : ofs_txt << setprecision(8) << rms_mtx2(0, 0) << endl;
3276 43 : ofs_txt << endl;
3277 43 : ofs_txt << "Number of clipped channels = "
3278 43 : << final_mask_after_clipping[ipol] - final_mask[ipol] << endl;
3279 43 : ofs_txt << endl;
3280 43 : ofs_txt << "------------------------------------------------------"
3281 43 : << endl;
3282 43 : ofs_txt << endl;
3283 : }
3284 : }
3285 :
3286 86 : if (write_baseline_csv) {
3287 12 : for (size_t ipol = 0; ipol < num_pol; ++ipol) {
3288 8 : if (apply_mtx2(ipol, 0) == false) continue;
3289 :
3290 6 : ofs_csv << (uInt)scans[irow] << ',' << (uInt)beams[irow] << ','
3291 6 : << (uInt)data_spw[irow] << ',' << ipol << ','
3292 6 : << setprecision(12) << times[irow] << ',';
3293 6 : ofs_csv << '[';
3294 13 : for (size_t imasklist = 0; imasklist < num_masklist_max / 2; ++imasklist) {
3295 7 : if (imasklist == 0) {
3296 6 : ofs_csv << '[' << masklist_mtx2(ipol, 2 * imasklist) << ';'
3297 6 : << masklist_mtx2(ipol, 2 * imasklist + 1) << ']';
3298 : }
3299 7 : if (imasklist >= 1
3300 8 : && (0 != masklist_mtx2(ipol, 2 * imasklist)
3301 1 : && 0 != masklist_mtx2(ipol, 2 * imasklist + 1))) {
3302 1 : ofs_csv << ";[" << masklist_mtx2(ipol, 2 * imasklist) << ';'
3303 1 : << masklist_mtx2(ipol, 2 * imasklist + 1) << ']';
3304 : }
3305 : }
3306 :
3307 6 : ofs_csv << ']' << ',';
3308 18 : Matrix<uInt> bltype_mtx2 = bltype_mtx[0][ipol];
3309 12 : string bltype_name;
3310 6 : if (bltype_mtx2(0, 0) == (uInt)0) {
3311 3 : bltype_name = "poly";
3312 3 : } else if (bltype_mtx2(0, 0) == (uInt)1) {
3313 2 : bltype_name = "chebyshev";
3314 1 : } else if (bltype_mtx2(0, 0) == (uInt)2) {
3315 1 : bltype_name = "cspline";
3316 : }
3317 :
3318 12 : Matrix<Int> fpar_mtx2 = fpar_mtx;
3319 12 : Matrix<Float> coeff_mtx2 = coeff_mtx;
3320 6 : ofs_csv << bltype_name.c_str() << ',' << fpar_mtx2(ipol, 0)
3321 6 : << ',';
3322 18 : for (size_t icoeff = 0; icoeff < numcoeff[ipol]; ++icoeff) {
3323 12 : ofs_csv << setprecision(8) << coeff_mtx2(ipol, icoeff) << ',';
3324 : }
3325 12 : Matrix<Float> rms_mtx2 = rms_mtx;
3326 6 : ofs_csv << setprecision(8) << rms_mtx2(ipol, 0) << ',';
3327 6 : ofs_csv << final_mask_after_clipping[ipol] - final_mask[ipol];
3328 6 : ofs_csv << endl;
3329 : }
3330 : }
3331 : } // end of chunk row loop
3332 : // write back data and flag cube to VisBuffer
3333 219 : if (update_weight) {
3334 6 : sdh_->fillCubeToOutputMs(vb, data_chunk, &flag_chunk, &weight_matrix);
3335 : } else {
3336 213 : sdh_->fillCubeToOutputMs(vb, data_chunk, &flag_chunk);
3337 : }
3338 219 : } // end of vi loop
3339 : } // end of chunk loop
3340 :
3341 81 : if (write_baseline_csv) {
3342 27 : ofs_csv.close();
3343 : }
3344 81 : if (write_baseline_text) {
3345 12 : ofs_txt.close();
3346 : }
3347 81 : if (write_baseline_table) {
3348 21 : bt->save(bloutputname_table);
3349 21 : delete bt;
3350 : }
3351 :
3352 81 : finalize_process();
3353 : // destroy baseline contexts
3354 81 : map<size_t const, std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> >::iterator ctxiter = context_reservoir.begin();
3355 289 : while (ctxiter != context_reservoir.end()) {
3356 208 : destroy_baseline_contexts (context_reservoir[(*ctxiter).first]);
3357 208 : ++ctxiter;
3358 : }
3359 81 : } //end subtractBaselineVariable
3360 :
3361 3504 : list<pair<size_t, size_t>> SingleDishMS::findLineAndGetRanges(size_t const num_data,
3362 : float const* data,
3363 : bool * mask,
3364 : float const threshold,
3365 : int const avg_limit,
3366 : int const minwidth,
3367 : vector<int> const& edge,
3368 : bool const invert) {
3369 : // input value check
3370 3504 : AlwaysAssert(minwidth > 0, AipsError);
3371 3504 : AlwaysAssert(avg_limit >= 0, AipsError);
3372 3504 : size_t max_iteration = 10;
3373 3504 : size_t maxwidth = num_data;
3374 3504 : AlwaysAssert(maxwidth > static_cast<size_t>(minwidth), AipsError);
3375 : // edge handling
3376 3504 : pair<size_t, size_t> lf_edge;
3377 3504 : if (edge.size() == 0) {
3378 0 : lf_edge = pair<size_t, size_t>(0, 0);
3379 3504 : } else if (edge.size() == 1) {
3380 2 : AlwaysAssert(edge[0] >= 0, AipsError);
3381 2 : lf_edge = pair<size_t, size_t>(static_cast<size_t>(edge[0]), static_cast<size_t>(edge[0]));
3382 : } else {
3383 3502 : AlwaysAssert(edge[0] >= 0 && edge[1] >= 0, AipsError);
3384 3502 : lf_edge = pair<size_t, size_t>(static_cast<size_t>(edge[0]), static_cast<size_t>(edge[1]));
3385 : }
3386 : // line detection
3387 : list<pair<size_t, size_t>> line_ranges = linefinder::MADLineFinder(num_data,
3388 : data, mask, threshold, max_iteration, static_cast<size_t>(minwidth),
3389 3504 : maxwidth, static_cast<size_t>(avg_limit), lf_edge);
3390 : // debug output
3391 10512 : LogIO os(_ORIGIN);
3392 3504 : os << LogIO::DEBUG1 << line_ranges.size() << " lines found: ";
3393 12762 : for (list<pair<size_t, size_t>>::iterator iter = line_ranges.begin();
3394 22020 : iter != line_ranges.end(); ++iter) {
3395 9258 : os << "[" << (*iter).first << ", " << (*iter).second << "] ";
3396 : }
3397 3504 : os << LogIO::POST;
3398 3504 : if (invert) { // eliminate edge channels from output mask
3399 3484 : if (lf_edge.first > 0)
3400 3 : line_ranges.push_front(pair<size_t, size_t>(0, lf_edge.first - 1));
3401 3484 : if (lf_edge.second > 0)
3402 3 : line_ranges.push_back(
3403 6 : pair<size_t, size_t>(num_data - lf_edge.second, num_data - 1));
3404 : }
3405 7008 : return line_ranges;
3406 : }
3407 :
3408 3484 : void SingleDishMS::findLineAndGetMask(size_t const num_data,
3409 : float const* data,
3410 : bool const* in_mask,
3411 : float const threshold,
3412 : int const avg_limit,
3413 : int const minwidth,
3414 : vector<int> const& edge,
3415 : bool const invert,
3416 : bool* out_mask) {
3417 : // copy input mask to output mask vector if necessary
3418 3484 : if (in_mask != out_mask) {
3419 0 : for (size_t i = 0; i < num_data; ++i) {
3420 0 : out_mask[i] = in_mask[i];
3421 : }
3422 : }
3423 : // line finding
3424 : list<pair<size_t, size_t>> line_ranges
3425 : = findLineAndGetRanges(num_data, data, out_mask, threshold,
3426 6968 : avg_limit, minwidth, edge, invert);
3427 : // line mask creation (do not initialize in case of baseline mask)
3428 3484 : linefinder::getMask(num_data, out_mask, line_ranges, invert, !invert);
3429 3484 : }
3430 :
3431 160 : void SingleDishMS::smooth(string const &kernelType,
3432 : float const kernelWidth,
3433 : string const &columnName,
3434 : string const &outMSName) {
3435 480 : LogIO os(_ORIGIN);
3436 : os << "Input parameter summary:" << endl << " kernelType = " << kernelType
3437 : << endl << " kernelWidth = " << kernelWidth << endl
3438 : << " columnName = " << columnName << endl << " outMSName = "
3439 160 : << outMSName << LogIO::POST;
3440 :
3441 : // Initialization
3442 160 : doSmoothing_ = true;
3443 160 : prepare_for_process(columnName, outMSName);
3444 :
3445 : // configure smoothing
3446 157 : sdh_->setSmoothing(kernelType, kernelWidth);
3447 157 : sdh_->initializeSmoothing();
3448 :
3449 : // get VI/VB2 access
3450 157 : vi::VisibilityIterator2 *visIter = sdh_->getVisIter();
3451 157 : visIter->setRowBlocking(kNRowBlocking);
3452 157 : vi::VisBuffer2 *vb = visIter->getVisBuffer();
3453 :
3454 157 : double startTime = gettimeofday_sec();
3455 :
3456 334 : for (visIter->originChunks(); visIter->moreChunks(); visIter->nextChunk()) {
3457 866 : for (visIter->origin(); visIter->more(); visIter->next()) {
3458 689 : sdh_->fillOutputMs(vb);
3459 : }
3460 : }
3461 :
3462 157 : double endTime = gettimeofday_sec();
3463 : os << LogIO::DEBUGGING
3464 : << "Elapsed time for VI/VB loop: " << endTime - startTime << " sec"
3465 157 : << LogIO::POST;
3466 :
3467 : // Finalization
3468 157 : finalize_process();
3469 157 : }
3470 :
3471 29 : void SingleDishMS::atmcor(Record const &config, string const &columnName, string const &outMSName) {
3472 87 : LogIO os(_ORIGIN);
3473 : os << LogIO::DEBUGGING
3474 : << "Input parameter summary:" << endl
3475 : << " columnName = " << columnName << endl << " outMSName = "
3476 29 : << outMSName << LogIO::POST;
3477 :
3478 : // Initialization
3479 29 : doAtmCor_ = true;
3480 29 : atmCorConfig_ = config;
3481 29 : os << LogIO::DEBUGGING << "config summry:";
3482 29 : atmCorConfig_.print(os.output(), 25, " ");
3483 29 : os << LogIO::POST;
3484 58 : Block<Int> sortCols(4);
3485 29 : sortCols[0] = MS::OBSERVATION_ID;
3486 29 : sortCols[1] = MS::ARRAY_ID;
3487 29 : sortCols[2] = MS::FEED1;
3488 29 : sortCols[3] = MS::DATA_DESC_ID;
3489 29 : prepare_for_process(columnName, outMSName, sortCols, False);
3490 :
3491 : // get VI/VB2 access
3492 27 : vi::VisibilityIterator2 *visIter = sdh_->getVisIter();
3493 : // for parallel processing: set row blocking (common multiple of 3 and 4)
3494 : // TODO: optimize row blocking size
3495 27 : constexpr rownr_t kNrowBlocking = 360u;
3496 81 : std::vector<Int> antenna1 = ScalarColumn<Int>(visIter->ms(), "ANTENNA1").getColumn().tovector();
3497 27 : std::sort(antenna1.begin(), antenna1.end());
3498 27 : auto const result = std::unique(antenna1.begin(), antenna1.end());
3499 27 : Int const nAntennas = std::distance(antenna1.begin(), result);
3500 27 : visIter->setRowBlocking(kNrowBlocking * nAntennas);
3501 : os << "There are " << nAntennas << " antennas in MAIN table. "
3502 : << "Set row-blocking size " << kNrowBlocking * nAntennas
3503 27 : << LogIO::POST;
3504 27 : vi::VisBuffer2 *vb = visIter->getVisBuffer();
3505 :
3506 27 : double startTime = gettimeofday_sec();
3507 :
3508 76 : for (visIter->originChunks(); visIter->moreChunks(); visIter->nextChunk()) {
3509 188 : for (visIter->origin(); visIter->more(); visIter->next()) {
3510 139 : sdh_->fillOutputMs(vb);
3511 : }
3512 : }
3513 :
3514 26 : double endTime = gettimeofday_sec();
3515 : os << LogIO::DEBUGGING
3516 : << "Elapsed time for VI/VB loop: " << endTime - startTime << " sec"
3517 26 : << LogIO::POST;
3518 :
3519 : // Finalization
3520 26 : finalize_process();
3521 26 : }
3522 :
3523 :
3524 5 : bool SingleDishMS::importAsap(string const &infile, string const &outfile, bool const parallel) {
3525 5 : bool status = true;
3526 : try {
3527 10 : SingleDishMSFiller<Scantable2MSReader> filler(infile, parallel);
3528 5 : filler.fill();
3529 5 : filler.save(outfile);
3530 0 : } catch (AipsError &e) {
3531 0 : LogIO os(_ORIGIN);
3532 0 : os << LogIO::SEVERE << "Exception occurred." << LogIO::POST;
3533 0 : os << LogIO::SEVERE << "Original Message: \n" << e.getMesg() << LogIO::POST;
3534 0 : os << LogIO::DEBUGGING << "Detailed Stack Trace: \n" << e.getStackTrace() << LogIO::POST;
3535 0 : status = false;
3536 0 : } catch (...) {
3537 0 : LogIO os(_ORIGIN);
3538 0 : os << LogIO::SEVERE << "Unknown exception occurred." << LogIO::POST;
3539 0 : status = false;
3540 : }
3541 5 : return status;
3542 : }
3543 :
3544 2 : bool SingleDishMS::importNRO(string const &infile, string const &outfile, bool const parallel) {
3545 2 : bool status = true;
3546 : try {
3547 4 : SingleDishMSFiller<NRO2MSReader> filler(infile, parallel);
3548 2 : filler.fill();
3549 2 : filler.save(outfile);
3550 0 : } catch (AipsError &e) {
3551 0 : LogIO os(_ORIGIN);
3552 0 : os << LogIO::SEVERE << "Exception occurred." << LogIO::POST;
3553 0 : os << LogIO::SEVERE << "Original Message: \n" << e.getMesg() << LogIO::POST;
3554 0 : os << LogIO::DEBUGGING << "Detailed Stack Trace: \n" << e.getStackTrace() << LogIO::POST;
3555 0 : status = false;
3556 0 : } catch (...) {
3557 0 : LogIO os(_ORIGIN);
3558 0 : os << LogIO::SEVERE << "Unknown exception occurred." << LogIO::POST;
3559 0 : status = false;
3560 : }
3561 2 : return status;
3562 : }
3563 :
3564 : } // End of casa namespace.
|
