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