#include <algorithm> #include <cstdlib> #include <iostream> #include <map> #include <string> #include <sys/time.h> #include <vector> #include <casacore/casa/Arrays/ArrayMath.h> #include <casacore/casa/Arrays/Vector.h> #include <casacore/casa/Containers/Allocator.h> #include <casacore/casa/Containers/Block.h> #include <casacore/casa/Logging/LogIO.h> #include <casacore/casa/Logging/LogOrigin.h> #include <casacore/casa/Quanta/MVTime.h> #include <casacore/casa/Utilities/Assert.h> #include <casacore/casa/Utilities/GenSort.h> #include <casacore/ms/MeasurementSets/MSSpectralWindow.h> #include <casacore/ms/MSSel/MSSelection.h> #include <casacore/ms/MSSel/MSSelectionTools.h> #include <msvis/MSVis/VisibilityIterator2.h> #include <msvis/MSVis/VisSetUtil.h> #include <casacore/scimath/Fitting/GenericL2Fit.h> #include <casacore/scimath/Fitting/NonLinearFitLM.h> #include <casacore/scimath/Functionals/CompiledFunction.h> #include <casacore/scimath/Functionals/CompoundFunction.h> #include <casacore/scimath/Functionals/Function.h> #include <casacore/scimath/Functionals/Gaussian1D.h> #include <casacore/scimath/Functionals/Lorentzian1D.h> #include <casacore/scimath/Mathematics/Convolver.h> #include <casacore/scimath/Mathematics/VectorKernel.h> #include <singledish/SingleDish/BaselineTable.h> #include <singledish/SingleDish/BLParameterParser.h> #include <singledish/SingleDish/LineFinder.h> #include <singledish/SingleDish/SingleDishMS.h> #include <stdcasa/StdCasa/CasacSupport.h> #include <casacore/tables/Tables/ScalarColumn.h> #include <casa_sakura/SakuraAlignedArray.h> // for importasap and importnro #include <singledishfiller/Filler/NRO2MSReader.h> #include <singledishfiller/Filler/Scantable2MSReader.h> #include <singledishfiller/Filler/SingleDishMSFiller.h> #define _ORIGIN LogOrigin("SingleDishMS", __func__, WHERE) namespace { // Max number of rows to get in each iteration constexpr casacore::Int kNRowBlocking = 1000; // Sinusoid constexpr int SinusoidWaveNumber_kUpperLimit = -999; // Weight constexpr size_t WeightIndex_kStddev = 0; constexpr size_t WeightIndex_kRms = 1; constexpr size_t WeightIndex_kNum = 2; double gettimeofday_sec() { struct timeval tv; gettimeofday(&tv, NULL); return tv.tv_sec + (double) tv.tv_usec * 1.0e-6; } using casa::vi::VisBuffer2; using casacore::Matrix; using casacore::Cube; using casacore::Float; using casacore::Complex; using casacore::AipsError; template<class CUBE_ACCESSOR> struct DataAccessorInterface { static void GetCube(VisBuffer2 const &vb, Cube<Float> &cube) { real(cube, CUBE_ACCESSOR::GetVisCube(vb)); } static void GetSlice(VisBuffer2 const &vb, size_t const iPol, Matrix<Float> &cubeSlice) { real(cubeSlice, CUBE_ACCESSOR::GetVisCube(vb).yzPlane(iPol)); } }; struct DataAccessor: public DataAccessorInterface<DataAccessor> { static Cube<Complex> GetVisCube(VisBuffer2 const &vb) { return vb.visCube(); } }; struct CorrectedDataAccessor: public DataAccessorInterface<CorrectedDataAccessor> { static Cube<Complex> GetVisCube(VisBuffer2 const &vb) { return vb.visCubeCorrected(); } }; struct FloatDataAccessor { static void GetCube(VisBuffer2 const &vb, Cube<Float> &cube) { cube = GetVisCube(vb); } static void GetSlice(VisBuffer2 const &vb, size_t const iPol, Matrix<Float> &cubeSlice) { cubeSlice = GetVisCube(vb).yzPlane(iPol); } private: static Cube<Float> GetVisCube(VisBuffer2 const &vb) { return vb.visCubeFloat(); } }; inline void GetCubeFromData(VisBuffer2 const &vb, Cube<Float> &cube) { DataAccessor::GetCube(vb, cube); } inline void GetCubeFromCorrected(VisBuffer2 const &vb, Cube<Float> &cube) { CorrectedDataAccessor::GetCube(vb, cube); } inline void GetCubeFromFloat(VisBuffer2 const &vb, Cube<Float> &cube) { FloatDataAccessor::GetCube(vb, cube); } inline void GetCubeDefault(VisBuffer2 const& /*vb*/, Cube<Float>& /*cube*/) { throw AipsError("Data accessor for VB2 is not properly configured."); } inline void compute_weight(size_t const num_data, float const data[/*num_data*/], bool const mask[/*num_data*/], std::vector<float>& weight) { for (size_t i = 0; i < WeightIndex_kNum; ++i) { weight[i] = 0.0; } int num_data_effective = 0; double sum = 0.0; double sum_sq = 0.0; for (size_t i = 0; i < num_data; ++i) { if (mask[i]) { num_data_effective++; sum += data[i]; sum_sq += data[i] * data[i]; } } if (num_data_effective > 0) { double factor = 1.0 / static_cast<double>(num_data_effective); double mean = sum * factor; double mean_sq = sum_sq * factor; std::vector<double> variance(WeightIndex_kNum); variance[WeightIndex_kStddev] = mean_sq - mean * mean; variance[WeightIndex_kRms] = mean_sq; auto do_compute_weight = [&](size_t idx) { if (variance[idx] > 0.0) { weight[idx] = static_cast<float>(1.0 / variance[idx]); } else { LogIO os(_ORIGIN); os << "Weight set to 0 for a bad data." << LogIO::WARN; } }; do_compute_weight(WeightIndex_kStddev); do_compute_weight(WeightIndex_kRms); } } } // anonymous namespace using namespace casacore; namespace casa { SingleDishMS::SingleDishMS() : msname_(""), sdh_(0) { initialize(); } SingleDishMS::SingleDishMS(string const& ms_name) : msname_(ms_name), sdh_(0) { LogIO os(_ORIGIN); initialize(); } SingleDishMS::~SingleDishMS() { if (sdh_) { delete sdh_; sdh_ = 0; } msname_ = ""; } void SingleDishMS::initialize() { in_column_ = MS::UNDEFINED_COLUMN; // out_column_ = MS::UNDEFINED_COLUMN; doSmoothing_ = false; doAtmCor_ = false; visCubeAccessor_ = GetCubeDefault; } bool SingleDishMS::close() { LogIO os(_ORIGIN); os << "Detaching from SingleDishMS" << LogIO::POST; if (sdh_) { delete sdh_; sdh_ = 0; } msname_ = ""; return true; } //////////////////////////////////////////////////////////////////////// ///// Common utility functions //////////////////////////////////////////////////////////////////////// void SingleDishMS::setSelection(Record const &selection, bool const verbose) { LogIO os(_ORIGIN); if (!selection_.empty()) // selection is set before os << "Discard old selection and setting new one." << LogIO::POST; if (selection.empty()) // empty selection is passed os << "Resetting selection." << LogIO::POST; selection_ = selection; // Verbose output bool any_selection(false); if (verbose && !selection_.empty()) { String timeExpr(""), antennaExpr(""), fieldExpr(""), spwExpr(""), uvDistExpr(""), taQLExpr(""), polnExpr(""), scanExpr(""), arrayExpr(""), obsExpr(""), intentExpr(""); timeExpr = get_field_as_casa_string(selection_, "timerange"); antennaExpr = get_field_as_casa_string(selection_, "antenna"); fieldExpr = get_field_as_casa_string(selection_, "field"); spwExpr = get_field_as_casa_string(selection_, "spw"); uvDistExpr = get_field_as_casa_string(selection_, "uvdist"); taQLExpr = get_field_as_casa_string(selection_, "taql"); polnExpr = get_field_as_casa_string(selection_, "correlation"); scanExpr = get_field_as_casa_string(selection_, "scan"); arrayExpr = get_field_as_casa_string(selection_, "array"); intentExpr = get_field_as_casa_string(selection_, "intent"); obsExpr = get_field_as_casa_string(selection_, "observation"); // selection Summary os << "[Selection Summary]" << LogIO::POST; if (obsExpr != "") { any_selection = true; os << "- Observation: " << obsExpr << LogIO::POST; } if (antennaExpr != "") { any_selection = true; os << "- Antenna: " << antennaExpr << LogIO::POST; } if (fieldExpr != "") { any_selection = true; os << "- Field: " << fieldExpr << LogIO::POST; } if (spwExpr != "") { any_selection = true; os << "- SPW: " << spwExpr << LogIO::POST; } if (polnExpr != "") { any_selection = true; os << "- Pol: " << polnExpr << LogIO::POST; } if (scanExpr != "") { any_selection = true; os << "- Scan: " << scanExpr << LogIO::POST; } if (timeExpr != "") { any_selection = true; os << "- Time: " << timeExpr << LogIO::POST; } if (intentExpr != "") { any_selection = true; os << "- Intent: " << intentExpr << LogIO::POST; } if (arrayExpr != "") { any_selection = true; os << "- Array: " << arrayExpr << LogIO::POST; } if (uvDistExpr != "") { any_selection = true; os << "- UVDist: " << uvDistExpr << LogIO::POST; } if (taQLExpr != "") { any_selection = true; os << "- TaQL: " << taQLExpr << LogIO::POST; } {// reindex Int ifield; ifield = selection_.fieldNumber(String("reindex")); if (ifield > -1) { Bool reindex = selection_.asBool(ifield); os << "- Reindex: " << (reindex ? "ON" : "OFF" ) << LogIO::POST; } } if (!any_selection) os << "No valid selection parameter is set." << LogIO::WARN; } } void SingleDishMS::setAverage(Record const &average, bool const verbose) { LogIO os(_ORIGIN); if (!average_.empty()) // average is set before os << "Discard old average and setting new one." << LogIO::POST; if (average.empty()) // empty average is passed os << "Resetting average." << LogIO::POST; average_ = average; if (verbose && !average_.empty()) { LogIO os(_ORIGIN); Int ifield; ifield = average_.fieldNumber(String("timeaverage")); os << "[Averaging Settings]" << LogIO::POST; if (ifield < 0 || !average_.asBool(ifield)) { os << "No averaging will be done" << LogIO::POST; return; } String timebinExpr(""), timespanExpr(""), tweightExpr(""); timebinExpr = get_field_as_casa_string(average_, "timebin"); timespanExpr = get_field_as_casa_string(average_, "timespan"); tweightExpr = get_field_as_casa_string(average_, "tweight"); if (timebinExpr != "") { os << "- Time bin: " << timebinExpr << LogIO::POST; } if (timespanExpr != "") { os << "- Time span: " << timespanExpr << LogIO::POST; } if (tweightExpr != "") { os << "- Averaging weight: " << tweightExpr << LogIO::POST; } } } void SingleDishMS::setPolAverage(Record const &average, bool const verbose) { LogIO os(_ORIGIN); if (!pol_average_.empty()) // polaverage is set before os << "Discard old average and setting new one." << LogIO::POST; if (average.empty()) // empty polaverage is passed os << "Resetting average." << LogIO::POST; pol_average_ = average; if (verbose && !pol_average_.empty()) { LogIO os(_ORIGIN); Int ifield; ifield = pol_average_.fieldNumber(String("polaverage")); os << "[Polarization Averaging Settings]" << LogIO::POST; if (ifield < 0 || !pol_average_.asBool(ifield)) { os << "No polarization averaging will be done" << LogIO::POST; return; } String polAverageModeExpr(""); polAverageModeExpr = get_field_as_casa_string(pol_average_, "polaveragemode"); if (polAverageModeExpr != "") { os << "- Mode: " << polAverageModeExpr << LogIO::POST; } } } String SingleDishMS::get_field_as_casa_string(Record const &in_data, string const &field_name) { Int ifield; ifield = in_data.fieldNumber(String(field_name)); if (ifield > -1) return in_data.asString(ifield); return ""; } bool SingleDishMS::prepare_for_process(string const &in_column_name, string const &out_ms_name) { // Sort by single dish default return prepare_for_process(in_column_name, out_ms_name, Block<Int>(), true); } bool SingleDishMS::prepare_for_process(string const &in_column_name, string const &out_ms_name, Block<Int> const &sortColumns, bool const addDefaultSortCols) { LogIO os(_ORIGIN); AlwaysAssert(msname_ != "", AipsError); // define a column to read data from string in_column_name_lower = in_column_name; std::transform( in_column_name_lower.begin(), in_column_name_lower.end(), in_column_name_lower.begin(), [](unsigned char c) {return std::tolower(c);} ); if (in_column_name_lower == "float_data") { in_column_ = MS::FLOAT_DATA; visCubeAccessor_ = GetCubeFromFloat; } else if (in_column_name_lower == "corrected") { in_column_ = MS::CORRECTED_DATA; visCubeAccessor_ = GetCubeFromCorrected; } else if (in_column_name_lower == "data") { in_column_ = MS::DATA; visCubeAccessor_ = GetCubeFromData; } else { throw(AipsError("Invalid data column name")); } // destroy SDMSManager if (sdh_) delete sdh_; // Configure record Record configure_param(selection_); format_selection(configure_param); configure_param.define("inputms", msname_); configure_param.define("outputms", out_ms_name); String in_name(in_column_name); in_name.upcase(); configure_param.define("datacolumn", in_name); // merge averaging parameters configure_param.merge(average_); // The other available keys // - buffermode, realmodelcol, usewtspectrum, tileshape, // - chanaverage, chanbin, useweights, // - ddistart, hanning // - regridms, phasecenter, restfreq, outframe, interpolation, nspw, // - mode, nchan, start, width, veltype, // - timeaverage, timebin, timespan, maxuvwdistance // smoothing configure_param.define("smoothFourier", doSmoothing_); // merge polarization averaging parameters configure_param.merge(pol_average_); // offline ATM correction configure_param.define("atmCor", doAtmCor_); configure_param.merge(atmCorConfig_); // Generate SDMSManager sdh_ = new SDMSManager(); // Configure SDMSManager sdh_->configure(configure_param); ostringstream oss; configure_param.print(oss); String str(oss.str()); os << LogIO::DEBUG1 << " Configuration Record " << LogIO::POST; os << LogIO::DEBUG1 << str << LogIO::POST; // Open the MS and select data sdh_->open(); sdh_->getOutputMs()->flush(); // set large timebin if not averaging Double timeBin; int exists = configure_param.fieldNumber("timebin"); if (exists < 0) { // Avoid TIME col being added to sort columns in MSIter::construct. // TIME is automatically added to sort column when // timebin is not 0, even if addDefaultSortCols=false. timeBin = 0.0; } else { String timebin_string; configure_param.get(exists, timebin_string); timeBin = casaQuantity(timebin_string).get("s").getValue(); Int ifield; ifield = configure_param.fieldNumber(String("timeaverage")); Bool average_time = ifield < 0 ? false : configure_param.asBool(ifield); if (timeBin < 0 || (average_time && timeBin == 0.0)) throw(AipsError("time bin should be positive")); } // set sort column sdh_->setSortColumns(sortColumns, addDefaultSortCols, timeBin); // Set up the Data Handler sdh_->setup(); return true; } void SingleDishMS::finalize_process() { initialize(); if (sdh_) { sdh_->close(); delete sdh_; sdh_ = 0; } } void SingleDishMS::format_selection(Record &selection) { // At this moment sdh_ is not supposed to be generated yet. LogIO os(_ORIGIN); // format spw String const spwSel(get_field_as_casa_string(selection, "spw")); selection.define("spw", spwSel == "" ? "*" : spwSel); // Select only auto-correlation String autoCorrSel(""); os << "Formatting antenna selection to select only auto-correlation" << LogIO::POST; String const antennaSel(get_field_as_casa_string(selection, "antenna")); os << LogIO::DEBUG1 << "Input antenna expression = " << antennaSel << LogIO::POST; if (antennaSel == "") { //Antenna selection is NOT set autoCorrSel = String("*&&&"); } else { //User defined antenna selection MeasurementSet MSobj = MeasurementSet(msname_); MeasurementSet* theMS = &MSobj; MSSelection theSelection; theSelection.setAntennaExpr(antennaSel); TableExprNode exprNode = theSelection.toTableExprNode(theMS); Vector<Int> ant1Vec = theSelection.getAntenna1List(); os << LogIO::DEBUG1 << ant1Vec.nelements() << " antenna(s) are selected. ID = "; for (uInt i = 0; i < ant1Vec.nelements(); ++i) { os << ant1Vec[i] << ", "; if (autoCorrSel != "") autoCorrSel += ";"; autoCorrSel += String::toString(ant1Vec[i]) + "&&&"; } os << LogIO::POST; } os << LogIO::DEBUG1 << "Auto-correlation selection string: " << autoCorrSel << LogIO::POST; selection.define("antenna", autoCorrSel); } void SingleDishMS::get_data_cube_float(vi::VisBuffer2 const &vb, Cube<Float> &data_cube) { // if (in_column_ == MS::FLOAT_DATA) { // data_cube = vb.visCubeFloat(); // } else { //need to convert Complex cube to Float // Cube<Complex> cdata_cube(data_cube.shape()); // if (in_column_ == MS::DATA) { // cdata_cube = vb.visCube(); // } else { //MS::CORRECTED_DATA // cdata_cube = vb.visCubeCorrected(); // } // // convert Complext to Float // convertArrayC2F(data_cube, cdata_cube); // } (*visCubeAccessor_)(vb, data_cube); } void SingleDishMS::convertArrayC2F(Array<Float> &to, Array<Complex> const &from) { if (to.nelements() == 0 && from.nelements() == 0) { return; } if (to.shape() != from.shape()) { throw(ArrayConformanceError("Array shape differs")); } Array<Complex>::const_iterator endFrom = from.end(); Array<Complex>::const_iterator iterFrom = from.begin(); for (Array<Float>::iterator iterTo = to.begin(); iterFrom != endFrom; ++iterFrom, ++iterTo) { *iterTo = iterFrom->real(); } } std::vector<string> SingleDishMS::split_string(string const &s, char delim) { std::vector<string> elems; string item; for (size_t i = 0; i < s.size(); ++i) { char ch = s.at(i); if (ch == delim) { if (!item.empty()) { elems.push_back(item); } item.clear(); } else { item += ch; } } if (!item.empty()) { elems.push_back(item); } return elems; } bool SingleDishMS::file_exists(string const &filename) { FILE *fp; if ((fp = fopen(filename.c_str(), "r")) == NULL) return false; fclose(fp); return true; } void SingleDishMS::parse_spw(string const &in_spw, Vector<Int> &rec_spw, Matrix<Int> &rec_chan, Vector<size_t> &nchan, Vector<Vector<Bool> > &mask, Vector<bool> &nchan_set) { Record selrec = sdh_->getSelRec(in_spw); rec_spw = selrec.asArrayInt("spw"); rec_chan = selrec.asArrayInt("channel"); nchan.resize(rec_spw.nelements()); mask.resize(rec_spw.nelements()); nchan_set.resize(rec_spw.nelements()); for (size_t i = 0; i < nchan_set.nelements(); ++i) { nchan_set(i) = false; } } void SingleDishMS::get_nchan_and_mask(Vector<Int> const &rec_spw, Vector<Int> const &data_spw, Matrix<Int> const &rec_chan, size_t const num_chan, Vector<size_t> &nchan, Vector<Vector<Bool> > &mask, Vector<bool> &nchan_set, bool &new_nchan) { new_nchan = false; for (size_t i = 0; i < rec_spw.nelements(); ++i) { //get nchan by spwid and set to nchan[] for (size_t j = 0; j < data_spw.nelements(); ++j) { if ((!nchan_set(i)) && (data_spw(j) == rec_spw(i))) { bool found = false; for (size_t k = 0; k < nchan.nelements(); ++k) { if (!nchan_set(k)) continue; if (nchan(k) == num_chan) found = true; } if (!found) { new_nchan = true; } nchan(i) = num_chan; nchan_set(i) = true; break; } } if (!nchan_set(i)) continue; mask(i).resize(nchan(i)); // generate mask get_mask_from_rec(rec_spw(i), rec_chan, mask(i), true); } } void SingleDishMS::get_mask_from_rec(Int spwid, Matrix<Int> const &rec_chan, Vector<Bool> &mask, bool initialize) { if (initialize) { for (size_t j = 0; j < mask.nelements(); ++j) { mask(j) = false; } } //construct a list of (start, end, stride, start, end, stride, ...) //from rec_chan for the spwid std::vector<uInt> edge; edge.clear(); for (size_t j = 0; j < rec_chan.nrow(); ++j) { if (rec_chan.row(j)(0) == spwid) { edge.push_back(rec_chan.row(j)(1)); edge.push_back(rec_chan.row(j)(2)); edge.push_back(rec_chan.row(j)(3)); } } //generate mask for (size_t j = 0; j < edge.size()-2; j += 3) { for (size_t k = edge[j]; k <= edge[j + 1] && k < mask.size(); k += edge[j + 2]) { mask(k) = true; } } } void SingleDishMS::get_masklist_from_mask(size_t const num_chan, bool const *mask, Vector<uInt> &masklist) { size_t const max_num_masklist = num_chan + 1; masklist.resize(max_num_masklist); // clear uInt last_idx = num_chan - 1; uInt num_masklist = 0; auto append = [&](uInt i){ masklist[num_masklist] = i; num_masklist++; }; if (mask[0]) { append(0); } for (uInt i = 1; i < last_idx; ++i) { if (!mask[i]) continue; // The following if-statements must be judged independently. // Don't put them together as a single statement by connecting with '||'. if (!mask[i - 1]) { append(i); } if (!mask[i + 1]) { append(i); } } if (mask[last_idx]) { if ((1 <= last_idx) && (!mask[last_idx - 1])) { append(last_idx); } append(last_idx); } masklist.resize(num_masklist, true); } void SingleDishMS::get_baseline_context(size_t const bltype, uint16_t order, size_t num_chan, Vector<size_t> const &nchan, Vector<bool> const &nchan_set, Vector<size_t> &ctx_indices, std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> &bl_contexts) { size_t idx = 0; bool found = false; for (size_t i = 0; i < nchan.nelements(); ++i) { if ((nchan_set[i])&&(nchan[i] == num_chan)) { idx = bl_contexts.size(); found = true; break; } } if (found) { for (size_t i = 0; i < nchan.nelements(); ++i) { if ((nchan_set[i])&&(nchan[i] == num_chan)) { ctx_indices[i] = idx; } } LIBSAKURA_SYMBOL(LSQFitContextFloat) *context; LIBSAKURA_SYMBOL(Status) status = LIBSAKURA_SYMBOL(Status_kNG); if ((bltype == BaselineType_kPolynomial)||(bltype == BaselineType_kChebyshev)) { status = LIBSAKURA_SYMBOL(CreateLSQFitContextPolynomialFloat)(static_cast<LIBSAKURA_SYMBOL(LSQFitType)>(bltype), static_cast<uint16_t>(order), num_chan, &context); } else if (bltype == BaselineType_kCubicSpline) { status = LIBSAKURA_SYMBOL(CreateLSQFitContextCubicSplineFloat)(static_cast<uint16_t>(order), num_chan, &context); } else if (bltype == BaselineType_kSinusoid) { status = LIBSAKURA_SYMBOL(CreateLSQFitContextSinusoidFloat)(static_cast<uint16_t>(order), num_chan, &context); } check_sakura_status("sakura_CreateLSQFitContextFloat", status); bl_contexts.push_back(context); } } void SingleDishMS::get_baseline_context(size_t const bltype, uint16_t order, size_t num_chan, size_t ispw, Vector<size_t> &ctx_indices, std::vector<size_t> & ctx_nchans, std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> &bl_contexts) { AlwaysAssert(bl_contexts.size() == ctx_nchans.size() || bl_contexts.size() == ctx_nchans.size()-1 , AipsError); size_t idx = 0; bool found = false; for (size_t i = 0; i < bl_contexts.size(); ++i) { if (ctx_nchans[i] == num_chan) { idx = i; found = true; break; } } if (found) { // contexts with the valid number of channels already exists. // just update idx to bl_contexts and return. ctx_indices[ispw] = idx; return; } // contexts with the number of channels is not yet in bl_contexts. // Need to create a new context. ctx_indices[ispw] = bl_contexts.size(); LIBSAKURA_SYMBOL(LSQFitContextFloat) *context; LIBSAKURA_SYMBOL(Status) status = LIBSAKURA_SYMBOL(Status_kNG); if ((bltype == BaselineType_kPolynomial)||(bltype == BaselineType_kChebyshev)) { status = LIBSAKURA_SYMBOL(CreateLSQFitContextPolynomialFloat)(static_cast<LIBSAKURA_SYMBOL(LSQFitType)>(bltype), static_cast<uint16_t>(order), num_chan, &context); } else if (bltype == BaselineType_kCubicSpline) { status = LIBSAKURA_SYMBOL(CreateLSQFitContextCubicSplineFloat)(static_cast<uint16_t>(order), num_chan, &context); } else if (bltype == BaselineType_kSinusoid) { status = LIBSAKURA_SYMBOL(CreateLSQFitContextSinusoidFloat)(static_cast<uint16_t>(order), num_chan, &context); } check_sakura_status("sakura_CreateLSQFitContextFloat", status); bl_contexts.push_back(context); if (ctx_nchans.size() != bl_contexts.size()) { ctx_nchans.push_back(num_chan); } AlwaysAssert(bl_contexts.size() == ctx_nchans.size(), AipsError); } void SingleDishMS::destroy_baseline_contexts(std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> &bl_contexts) { LIBSAKURA_SYMBOL(Status) status; for (size_t i = 0; i < bl_contexts.size(); ++i) { status = LIBSAKURA_SYMBOL(DestroyLSQFitContextFloat)(bl_contexts[i]); check_sakura_status("sakura_DestoyBaselineContextFloat", status); } } void SingleDishMS::check_sakura_status(string const &name, LIBSAKURA_SYMBOL(Status) const status) { if (status == LIBSAKURA_SYMBOL(Status_kOK)) return; ostringstream oss; oss << name << "() failed -- "; if (status == LIBSAKURA_SYMBOL(Status_kNG)) { oss << "NG"; } else if (status == LIBSAKURA_SYMBOL(Status_kInvalidArgument)) { oss << "InvalidArgument"; } else if (status == LIBSAKURA_SYMBOL(Status_kNoMemory)) { oss << "NoMemory"; } else if (status == LIBSAKURA_SYMBOL(Status_kUnknownError)) { oss << "UnknownError"; } throw(AipsError(oss.str())); } void SingleDishMS::check_baseline_status(LIBSAKURA_SYMBOL(LSQFitStatus) const bl_status) { if (bl_status != LIBSAKURA_SYMBOL(LSQFitStatus_kOK)) { throw(AipsError("baseline fitting isn't successful.")); } } void SingleDishMS::get_spectrum_from_cube(Cube<Float> &data_cube, size_t const row, size_t const plane, size_t const num_data, float out_data[]) { AlwaysAssert(static_cast<size_t>(data_cube.ncolumn()) == num_data, AipsError); for (size_t i = 0; i < num_data; ++i) out_data[i] = static_cast<float>(data_cube(plane, i, row)); } void SingleDishMS::set_spectrum_to_cube(Cube<Float> &data_cube, size_t const row, size_t const plane, size_t const num_data, float *in_data) { AlwaysAssert(static_cast<size_t>(data_cube.ncolumn()) == num_data, AipsError); for (size_t i = 0; i < num_data; ++i) data_cube(plane, i, row) = static_cast<Float>(in_data[i]); } void SingleDishMS::get_weight_matrix(vi::VisBuffer2 const &vb, Matrix<Float> &weight_matrix) { weight_matrix = vb.weight(); } void SingleDishMS::set_weight_to_matrix(Matrix<Float> &weight_matrix, size_t const row, size_t const plane, float in_weight) { weight_matrix(plane, row) = static_cast<Float>(in_weight); } void SingleDishMS::get_flag_cube(vi::VisBuffer2 const &vb, Cube<Bool> &flag_cube) { flag_cube = vb.flagCube(); } void SingleDishMS::get_flag_from_cube(Cube<Bool> &flag_cube, size_t const row, size_t const plane, size_t const num_flag, bool out_flag[]) { AlwaysAssert(static_cast<size_t>(flag_cube.ncolumn()) == num_flag, AipsError); for (size_t i = 0; i < num_flag; ++i) out_flag[i] = static_cast<bool>(flag_cube(plane, i, row)); } void SingleDishMS::set_flag_to_cube(Cube<Bool> &flag_cube, size_t const row, size_t const plane, size_t const num_flag, bool *in_flag) { AlwaysAssert(static_cast<size_t>(flag_cube.ncolumn()) == num_flag, AipsError); for (size_t i = 0; i < num_flag; ++i) flag_cube(plane, i, row) = static_cast<Bool>(in_flag[i]); } void SingleDishMS::flag_spectrum_in_cube(Cube<Bool> &flag_cube, size_t const row, size_t const plane) { uInt const num_flag = flag_cube.ncolumn(); for (uInt ichan = 0; ichan < num_flag; ++ichan) flag_cube(plane, ichan, row) = true; } bool SingleDishMS::allchannels_flagged(size_t const num_flag, bool const* flag) { bool res = true; for (size_t i = 0; i < num_flag; ++i) { if (!flag[i]) { res = false; break; } } return res; } size_t SingleDishMS::NValidMask(size_t const num_mask, bool const* mask) { std::size_t nvalid = 0; // the assertion lines had better be replaced with static_assert when c++11 is supported AlwaysAssert(static_cast<std::size_t>(true) == 1, AipsError); AlwaysAssert(static_cast<std::size_t>(false) == 0, AipsError); for (size_t i = 0; i < num_mask; ++i) { nvalid += static_cast<std::size_t>(mask[i]); } return nvalid; } void SingleDishMS::split_bloutputname(string str) { char key = ','; vector<size_t> v; for (size_t i = 0; i < str.size(); ++i) { char target = str[i]; if (key == target) { v.push_back(i); } } //cout << "comma " << v.size() << endl; //cout << "v[1] " << v[1] << endl; //cout << "v.size()-1 " << v.size()-1 << endl; //cout << "v[1]+1 " << v[1]+1 << endl; //cout << "str.size()-v[1]-1 " << str.size()-v[1]-1 << endl; //cout << "str.substr(v[1]+1, str.size()-v[1]-1) " << str.substr(v[1]+1, str.size()-v[1]-1) << endl; string ss; bloutputname_csv.clear(); bloutputname_text.clear(); bloutputname_table.clear(); if (0 != v[0]) { bloutputname_csv = str.substr(0, v[0]); ss = str.substr(0, v[0]); } if (v[0] + 1 != v[1]) { bloutputname_text = str.substr(v[0] + 1, v[1] - v[0] - 1); } if (v[1] != str.size() - 1) { bloutputname_table = str.substr(v[1] + 1, str.size() - v[1] - 1); } } size_t SingleDishMS::get_num_coeff_bloutput(size_t const bltype, size_t order, size_t &num_coeff_max) { size_t num_coeff = 0; switch (bltype) { case BaselineType_kPolynomial: case BaselineType_kChebyshev: break; case BaselineType_kCubicSpline: num_coeff = order + 1; break; case BaselineType_kSinusoid: break; default: throw(AipsError("Unsupported baseline type.")); } if (num_coeff_max < num_coeff) { num_coeff_max = num_coeff; } return num_coeff; } vector<int> SingleDishMS::string_to_list(string const &wn_str, char const delim) { vector<int> wn_list; wn_list.clear(); vector<size_t> delim_position; delim_position.clear(); for (size_t i = 0; i < wn_str.size(); ++i) { if (wn_str[i] == delim) { delim_position.push_back(i); } } delim_position.push_back(wn_str.size()); size_t start_position = 0; for (size_t i = 0; i < delim_position.size(); ++i) { size_t end_position = delim_position[i]; size_t length = end_position - start_position; if (length > 0) { wn_list.push_back(std::atoi(wn_str.substr(start_position, length).c_str())); } start_position = end_position + 1; } return wn_list; } void SingleDishMS::get_effective_nwave(std::vector<int> const &addwn, std::vector<int> const &rejwn, int const wn_ulimit, std::vector<int> &effwn) { effwn.clear(); if (addwn.size() < 1) { throw AipsError("addwn has no elements."); } auto up_to_nyquist_limit = [&](std::vector<int> const &v){ return ((v.size() == 2) && (v[1] == SinusoidWaveNumber_kUpperLimit)); }; auto check_rejwn_add = [&](int const v){ bool add = (0 <= v) && (v <= wn_ulimit); // check v in range for (size_t i = 0; i < rejwn.size(); ++i) { if (rejwn[i] == v) { add = false; break; } } if (add) { effwn.push_back(v); } }; if (up_to_nyquist_limit(addwn)) { if (up_to_nyquist_limit(rejwn)) { if (addwn[0] < rejwn[0]) { for (int wn = addwn[0]; wn < rejwn[0]; ++wn) { if ((0 <= wn) && (wn <= wn_ulimit)) { effwn.push_back(wn); } } } else { throw(AipsError("No effective wave number given for sinusoidal fitting.")); } } else { for (int wn = addwn[0]; wn <= wn_ulimit; ++wn) { check_rejwn_add(wn); } } } else { if (up_to_nyquist_limit(rejwn)) { int maxwn = rejwn[0] - 1; if (maxwn < 0) { throw(AipsError("No effective wave number given for sinusoidal fitting.")); } for (size_t i = 0; i < addwn.size(); ++i) { if ((0 <= addwn[i]) && (addwn[i] <= maxwn)) { effwn.push_back(addwn[i]); } } } else { for (size_t i = 0; i < addwn.size(); ++i) { check_rejwn_add(addwn[i]); } } } if (effwn.size() == 0) { throw(AipsError("No effective wave number given for sinusoidal fitting.")); } } void SingleDishMS::finalise_effective_nwave(std::vector<int> const &blparam_eff_base, std::vector<int> const &blparam_exclude, int const &blparam_upperlimit, size_t const &num_chan, float const *spec, bool const *mask, bool const &applyfft, string const &fftmethod, string const &fftthresh_str, std::vector<size_t> &blparam_eff) { //why? blparam_eff.resize(blparam_eff_base.size()); copy(blparam_eff_base.begin(), blparam_eff_base.end(), blparam_eff.begin()); if (applyfft) { string fftthresh_attr; float fftthresh_sigma; int fftthresh_top; parse_fftthresh(fftthresh_str, fftthresh_attr, fftthresh_sigma, fftthresh_top); std::vector<int> blparam_fft; select_wavenumbers_via_fft(num_chan, spec, mask, fftmethod, fftthresh_attr, fftthresh_sigma, fftthresh_top, blparam_upperlimit, blparam_fft); merge_wavenumbers(blparam_eff_base, blparam_fft, blparam_exclude, blparam_eff); } } void SingleDishMS::parse_fftthresh(string const& fftthresh_str, string& fftthresh_attr, float& fftthresh_sigma, int& fftthresh_top) { size_t idx_sigma = fftthresh_str.find("sigma"); size_t idx_top = fftthresh_str.find("top"); if (idx_top == 0) { std::istringstream is(fftthresh_str.substr(3)); is >> fftthresh_top; fftthresh_attr = "top"; } else if (idx_sigma == fftthresh_str.size() - 5) { std::istringstream is(fftthresh_str.substr(0, fftthresh_str.size() - 5)); is >> fftthresh_sigma; fftthresh_attr = "sigma"; } else { bool is_number = true; for (size_t i = 0; i < fftthresh_str.size()-1; ++i) { char ch = (fftthresh_str.substr(i, 1).c_str())[0]; if (!(isdigit(ch) || (fftthresh_str.substr(i, 1) == "."))) { is_number = false; break; } } if (is_number) { std::istringstream is(fftthresh_str); is >> fftthresh_sigma; fftthresh_attr = "sigma"; } else { throw(AipsError("fftthresh has a wrong value")); } } } void SingleDishMS::select_wavenumbers_via_fft(size_t const num_chan, float const *spec, bool const *mask, string const &fftmethod, string const &fftthresh_attr, float const fftthresh_sigma, int const fftthresh_top, int const blparam_upperlimit, std::vector<int> &blparam_fft) { blparam_fft.clear(); std::vector<float> fourier_spec; if (fftmethod == "fft") { exec_fft(num_chan, spec, mask, false, true, fourier_spec); } else { throw AipsError("fftmethod must be 'fft' for now."); } // Anything except fft and 3.0 is not used? // top, sigma are not documented int fourier_spec_size = static_cast<int>(fourier_spec.size()); if (fftthresh_attr == "sigma") { float mean = 0.0; float mean2 = 0.0; for (int i = 0; i < fourier_spec_size; ++i) { mean += fourier_spec[i]; mean2 += fourier_spec[i] * fourier_spec[i]; } mean /= static_cast<float>(fourier_spec_size); mean2 /= static_cast<float>(fourier_spec_size); float thres = mean + fftthresh_sigma * float(sqrt(mean2 - mean * mean)); for (int i = 0; i < fourier_spec_size; ++i) { if ((i <= blparam_upperlimit)&&(thres <= fourier_spec[i])) { blparam_fft.push_back(i); } } } else if (fftthresh_attr == "top") { int i = 0; while (i < fftthresh_top) { float max = 0.0; int max_idx = 0; for (int j = 0; j < fourier_spec_size; ++j) { if (max < fourier_spec[j]) { max = fourier_spec[j]; max_idx = j; } } fourier_spec[max_idx] = 0.0; if (max_idx <= blparam_upperlimit) { blparam_fft.push_back(max_idx); ++i; } } } else { throw AipsError("fftthresh is wrong."); } } void SingleDishMS::exec_fft(size_t const num_chan, float const *in_spec, bool const *in_mask, bool const get_real_imag, bool const get_ampl_only, std::vector<float> &fourier_spec) { Vector<Float> spec; interpolate_constant(static_cast<int>(num_chan), in_spec, in_mask, spec); FFTServer<Float, Complex> ffts; Vector<Complex> fftres; ffts.fft0(fftres, spec); float norm = static_cast<float>(2.0/static_cast<double>(num_chan)); fourier_spec.clear(); if (get_real_imag) { for (size_t i = 0; i < fftres.size(); ++i) { fourier_spec.push_back(real(fftres[i]) * norm); fourier_spec.push_back(imag(fftres[i]) * norm); } } else { //not used? for (size_t i = 0; i < fftres.size(); ++i) { fourier_spec.push_back(abs(fftres[i]) * norm); if (!get_ampl_only) fourier_spec.push_back(arg(fftres[i])); } } } void SingleDishMS::interpolate_constant(int const num_chan, float const *in_spec, bool const *in_mask, Vector<Float> &spec) { spec.resize(num_chan); for (int i = 0; i < num_chan; ++i) { spec[i] = in_spec[i]; } int idx_left = -1; int idx_right = -1; bool masked_region = false; for (int i = 0; i < num_chan; ++i) { if (!in_mask[i]) { masked_region = true; idx_left = i; while (i < num_chan) { if (in_mask[i]) break; idx_right = i; ++i; } } if (masked_region) { // execute interpolation as the following criteria: // (1) for a masked region inside the spectrum, replace the spectral // values with the mean of those at the two channels just outside // the both edges of the masked region. // (2) for a masked region at the spectral edge, replace the values // with the one at the nearest non-masked channel. // (ZOH, but bilateral) Float interp = 0.0; int idx_left_next = idx_left - 1; int idx_right_next = idx_right + 1; if (idx_left_next < 0) { if (idx_right_next < num_chan) { interp = in_spec[idx_right_next]; } else { throw AipsError("Bad data. all channels are masked."); } } else { interp = in_spec[idx_left_next]; if (idx_right_next < num_chan) { interp = (interp + in_spec[idx_right_next]) / 2.0; } } if ((0 <= idx_left) && (idx_left < num_chan) && (0 <= idx_right) && (idx_right < num_chan)) { for (int j = idx_left; j <= idx_right; ++j) { spec[j] = interp; } } masked_region = false; } } } void SingleDishMS::merge_wavenumbers(std::vector<int> const &blparam_eff_base, std::vector<int> const &blparam_fft, std::vector<int> const &blparam_exclude, std::vector<size_t> &blparam_eff) { for (size_t i = 0; i < blparam_fft.size(); ++i) { bool found = false; for (size_t j = 0; j < blparam_eff_base.size(); ++j) { if (blparam_eff_base[j] == blparam_fft[i]) { found = true; break; } } if (!found) { //new value to add //but still need to check if it is to be excluded bool found_exclude = false; for (size_t j = 0; j < blparam_exclude.size(); ++j) { if (blparam_exclude[j] == blparam_fft[i]) { found_exclude = true; break; } } if (!found_exclude) { blparam_eff.push_back(blparam_fft[i]); } } } if (1 < blparam_eff.size()) { sort(blparam_eff.begin(), blparam_eff.end()); unique(blparam_eff.begin(), blparam_eff.end()); } } template<typename Func0, typename Func1, typename Func2, typename Func3> void SingleDishMS::doSubtractBaseline(string const& in_column_name, string const& out_ms_name, string const& out_bloutput_name, bool const& do_subtract, string const& in_spw, bool const& update_weight, string const& sigma_value, LIBSAKURA_SYMBOL(Status)& status, std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> &bl_contexts, size_t const bltype, vector<int> const& blparam, vector<int> const& blparam_exclude, bool const& applyfft, string const& fftmethod, string const& fftthresh, float const clip_threshold_sigma, int const num_fitting_max, bool const linefinding, float const threshold, int const avg_limit, int const minwidth, vector<int> const& edge, Func0 func0, Func1 func1, Func2 func2, Func3 func3, LogIO os) { os << LogIO::DEBUG2 << "Calling SingleDishMS::doSubtractBaseline " << LogIO::POST; // in_ms = out_ms // in_column = [FLOAT_DATA|DATA|CORRECTED_DATA], out_column=new MS // no iteration is necessary for the processing. // procedure // 1. iterate over MS // 2. pick single spectrum from in_column (this is not actually necessary for simple scaling but for exibision purpose) // 3. fit baseline to each spectrum and subtract it // 4. put single spectrum (or a block of spectra) to out_column Block<Int> columns(1); columns[0] = MS::DATA_DESC_ID; prepare_for_process(in_column_name, out_ms_name, columns, false); vi::VisibilityIterator2 *vi = sdh_->getVisIter(); vi->setRowBlocking(kNRowBlocking); vi::VisBuffer2 *vb = vi->getVisBuffer(); split_bloutputname(out_bloutput_name); bool write_baseline_csv = (bloutputname_csv != ""); bool write_baseline_text = (bloutputname_text != ""); bool write_baseline_table = (bloutputname_table != ""); ofstream ofs_csv; ofstream ofs_txt; BaselineTable *bt = 0; if (write_baseline_csv) { ofs_csv.open(bloutputname_csv.c_str()); } if (write_baseline_text) { ofs_txt.open(bloutputname_text.c_str(), std::ios_base::out | std::ios_base::app); } if (write_baseline_table) { bt = new BaselineTable(vi->ms()); } Vector<Int> recspw; Matrix<Int> recchan; Vector<size_t> nchan; Vector<Vector<Bool> > in_mask; Vector<bool> nchan_set; parse_spw(in_spw, recspw, recchan, nchan, in_mask, nchan_set); Vector<size_t> ctx_indices(nchan.nelements(), 0ul); std::vector<int> blparam_eff_base; auto wn_ulimit_by_rejwn = [&](){ return ((blparam_exclude.size() == 2) && (blparam_exclude[1] == SinusoidWaveNumber_kUpperLimit)); }; std::vector<float> weight(WeightIndex_kNum); size_t const var_index = (sigma_value == "stddev") ? WeightIndex_kStddev : WeightIndex_kRms; for (vi->originChunks(); vi->moreChunks(); vi->nextChunk()) { for (vi->origin(); vi->more(); vi->next()) { Vector<Int> scans = vb->scan(); Vector<Double> times = vb->time(); Vector<Int> beams = vb->feed1(); Vector<Int> antennas = vb->antenna1(); Vector<Int> data_spw = vb->spectralWindows(); size_t const num_chan = static_cast<size_t>(vb->nChannels()); size_t const num_pol = static_cast<size_t>(vb->nCorrelations()); size_t const num_row = static_cast<size_t>(vb->nRows()); Cube<Float> data_chunk(num_pol, num_chan, num_row, Array<Float>::uninitialized); SakuraAlignedArray<float> spec(num_chan); Cube<Bool> flag_chunk(num_pol, num_chan, num_row, Array<Bool>::uninitialized); SakuraAlignedArray<bool> flag(num_chan); SakuraAlignedArray<bool> mask(num_chan); SakuraAlignedArray<bool> mask_after_clipping(num_chan); float *spec_data = spec.data(); bool *flag_data = flag.data(); bool *mask_data = mask.data(); bool *mask_after_clipping_data = mask_after_clipping.data(); Matrix<Float> weight_matrix(num_pol, num_row, Array<Float>::uninitialized); auto get_wavenumber_upperlimit = [&](){ return static_cast<int>(num_chan) / 2 - 1; }; uInt final_mask[num_pol]; uInt final_mask_after_clipping[num_pol]; final_mask[0] = 0; final_mask[1] = 0; final_mask_after_clipping[0] = 0; final_mask_after_clipping[1] = 0; bool new_nchan = false; get_nchan_and_mask(recspw, data_spw, recchan, num_chan, nchan, in_mask, nchan_set, new_nchan); if (new_nchan) { int blparam_max = blparam[blparam.size() - 1]; if (bltype == BaselineType_kSinusoid) { int nwave_ulimit = get_wavenumber_upperlimit(); get_effective_nwave(blparam, blparam_exclude, nwave_ulimit, blparam_eff_base); blparam_max = blparam_eff_base[blparam_eff_base.size() - 1]; } if ((bltype != BaselineType_kSinusoid) || (!applyfft) || wn_ulimit_by_rejwn()) { get_baseline_context(bltype, static_cast<uint16_t>(blparam_max), num_chan, nchan, nchan_set, ctx_indices, bl_contexts); } } else { int last_nchan_set_idx = nchan_set.size() - 1; for (int i = nchan_set.size()-1; i >= 0; --i) { if (nchan_set[i]) break; --last_nchan_set_idx; } if (0 < last_nchan_set_idx) { for (int i = 0; i < last_nchan_set_idx; ++i) { if (nchan[i] == nchan[last_nchan_set_idx]) { ctx_indices[last_nchan_set_idx] = ctx_indices[i]; break; } } } } // get data/flag cubes (npol*nchan*nrow) from VisBuffer get_data_cube_float(*vb, data_chunk); get_flag_cube(*vb, flag_chunk); // get weight matrix (npol*nrow) from VisBuffer if (update_weight) { get_weight_matrix(*vb, weight_matrix); } // loop over MS rows for (size_t irow = 0; irow < num_row; ++irow) { size_t idx = 0; for (size_t ispw = 0; ispw < recspw.nelements(); ++ispw) { if (data_spw[irow] == recspw[ispw]) { idx = ispw; break; } } //prepare variables for writing baseline table Array<Bool> apply_mtx(IPosition(2, num_pol, 1), true); Array<uInt> bltype_mtx(IPosition(2, num_pol, 1), (uInt)bltype); //Array<Int> fpar_mtx(IPosition(2, num_pol, 1), (Int)blparam[blparam.size()-1]); std::vector<std::vector<size_t> > fpar_mtx_tmp(num_pol); std::vector<std::vector<double> > ffpar_mtx_tmp(num_pol); std::vector<std::vector<uInt> > masklist_mtx_tmp(num_pol); std::vector<std::vector<double> > coeff_mtx_tmp(num_pol); Array<Float> rms_mtx(IPosition(2, num_pol, 1), (Float)0); Array<Float> cthres_mtx(IPosition(2, num_pol, 1), Array<Float>::uninitialized); Array<uInt> citer_mtx(IPosition(2, num_pol, 1), Array<uInt>::uninitialized); Array<Bool> uself_mtx(IPosition(2, num_pol, 1), Array<Bool>::uninitialized); Array<Float> lfthres_mtx(IPosition(2, num_pol, 1), Array<Float>::uninitialized); Array<uInt> lfavg_mtx(IPosition(2, num_pol, 1), Array<uInt>::uninitialized); Array<uInt> lfedge_mtx(IPosition(2, num_pol, 2), Array<uInt>::uninitialized); size_t num_apply_true = 0; size_t num_fpar_max = 0; size_t num_ffpar_max = 0; size_t num_masklist_max = 0; size_t num_coeff_max = 0; // loop over polarization for (size_t ipol = 0; ipol < num_pol; ++ipol) { // get a channel mask from data cube // (note that the variable 'mask' is flag in the next line // actually, then it will be converted to real mask when // taking AND with user-given mask info. this is just for // saving memory usage...) get_flag_from_cube(flag_chunk, irow, ipol, num_chan, flag_data); // skip spectrum if all channels flagged if (allchannels_flagged(num_chan, flag_data)) { apply_mtx[0][ipol] = false; continue; } // convert flag to mask by taking logical NOT of flag // and then operate logical AND with in_mask for (size_t ichan = 0; ichan < num_chan; ++ichan) { mask_data[ichan] = in_mask[idx][ichan] && (!(flag_data[ichan])); } // get a spectrum from data cube get_spectrum_from_cube(data_chunk, irow, ipol, num_chan, spec_data); // line finding. get baseline mask (invert=true) if (linefinding) { findLineAndGetMask(num_chan, spec_data, mask_data, threshold, avg_limit, minwidth, edge, true, mask_data); } std::vector<size_t> blparam_eff; size_t num_coeff; if (bltype == BaselineType_kSinusoid) { int nwave_ulimit = get_wavenumber_upperlimit(); finalise_effective_nwave(blparam_eff_base, blparam_exclude, nwave_ulimit, num_chan, spec_data, mask_data, applyfft, fftmethod, fftthresh, blparam_eff); size_t blparam_eff_size = blparam_eff.size(); if (blparam_eff[0] == 0) { num_coeff = blparam_eff_size * 2 - 1; } else { num_coeff = blparam_eff_size * 2; } } else if (bltype == BaselineType_kCubicSpline) { blparam_eff.resize(1); blparam_eff[0] = blparam[blparam.size() - 1]; num_coeff = blparam_eff[0] * 4; } else { // poly, chebyshev blparam_eff.resize(1); blparam_eff[0] = blparam[blparam.size() - 1]; status = LIBSAKURA_SYMBOL(GetNumberOfCoefficientsFloat)(bl_contexts[ctx_indices[idx]], blparam_eff[0], &num_coeff); check_sakura_status("sakura_GetNumberOfCoefficients", status); } // Final check of the valid number of channels size_t num_min = (bltype == BaselineType_kCubicSpline) ? blparam[blparam.size()-1] + 3 : num_coeff; if (NValidMask(num_chan, mask_data) < num_min) { flag_spectrum_in_cube(flag_chunk, irow, ipol); apply_mtx[0][ipol] = false; os << LogIO::WARN << "Too few valid channels to fit. Skipping Antenna " << antennas[irow] << ", Beam " << beams[irow] << ", SPW " << data_spw[irow] << ", Pol " << ipol << ", Time " << MVTime(times[irow] / 24. / 3600.).string(MVTime::YMD, 8) << LogIO::POST; continue; } // actual execution of single spectrum float rms; if (write_baseline_text || write_baseline_csv || write_baseline_table) { num_apply_true++; if (num_coeff_max < num_coeff) { num_coeff_max = num_coeff; } SakuraAlignedArray<double> coeff(num_coeff); double *coeff_data = coeff.data(); //---GetBestFitBaselineCoefficientsFloat()... //func0(ctx_indices[idx], num_chan, blparam_eff, spec_data, mask_data, num_coeff, coeff_data, mask_after_clipping_data, &rms); LIBSAKURA_SYMBOL(LSQFitContextFloat) *context = nullptr; if ((bltype != BaselineType_kSinusoid) || (!applyfft) || wn_ulimit_by_rejwn()) { context = bl_contexts[ctx_indices[idx]]; } func0(context, num_chan, blparam_eff, spec_data, mask_data, num_coeff, coeff_data, mask_after_clipping_data, &rms); for (size_t i = 0; i < num_chan; ++i) { if (mask_data[i] == false) { final_mask[ipol] += 1; } if (mask_after_clipping_data[i] == false) { final_mask_after_clipping[ipol] += 1; } } //set_array_for_bltable(fpar_mtx_tmp) size_t num_fpar = blparam_eff.size(); fpar_mtx_tmp[ipol].resize(num_fpar); if (num_fpar_max < num_fpar) { num_fpar_max = num_fpar; } fpar_mtx_tmp[ipol].resize(num_fpar); for (size_t ifpar = 0; ifpar < num_fpar; ++ifpar) { fpar_mtx_tmp[ipol][ifpar] = blparam_eff[ifpar]; } //---set_array_for_bltable(ffpar_mtx_tmp) func1(ipol, ffpar_mtx_tmp, num_ffpar_max); //set_array_for_bltable<double, Float>(ipol, num_coeff, coeff_data, coeff_mtx); coeff_mtx_tmp[ipol].resize(num_coeff); for (size_t icoeff = 0; icoeff < num_coeff; ++icoeff) { coeff_mtx_tmp[ipol][icoeff] = coeff_data[icoeff]; } Vector<uInt> masklist; get_masklist_from_mask(num_chan, mask_after_clipping_data, masklist); if (masklist.size() > num_masklist_max) { num_masklist_max = masklist.size(); } masklist_mtx_tmp[ipol].clear(); for (size_t imask = 0; imask < masklist.size(); ++imask) { masklist_mtx_tmp[ipol].push_back(masklist[imask]); } //---SubtractBaselineUsingCoefficientsFloat()... //func2(ctx_indices[idx], num_chan, fpar_mtx_tmp[ipol], spec_data, num_coeff, coeff_data); func2(context, num_chan, fpar_mtx_tmp[ipol], spec_data, num_coeff, coeff_data); rms_mtx[0][ipol] = rms; cthres_mtx[0][ipol] = clip_threshold_sigma; citer_mtx[0][ipol] = (uInt)num_fitting_max - 1; uself_mtx[0][ipol] = false; lfthres_mtx[0][ipol] = 0.0; lfavg_mtx[0][ipol] = 0; for (size_t iedge = 0; iedge < 2; ++iedge) { lfedge_mtx[iedge][ipol] = 0; } } else { //---SubtractBaselineFloat()... //func3(ctx_indices[idx], num_chan, blparam_eff, num_coeff, spec_data, mask_data, &rms); LIBSAKURA_SYMBOL(LSQFitContextFloat) *context = nullptr; if ((bltype != BaselineType_kSinusoid) || (!applyfft) || wn_ulimit_by_rejwn()) { context = bl_contexts[ctx_indices[idx]]; } func3(context, num_chan, blparam_eff, num_coeff, spec_data, mask_data, mask_after_clipping_data, &rms); } // set back a spectrum to data cube if (do_subtract) { set_spectrum_to_cube(data_chunk, irow, ipol, num_chan, spec_data); } if (update_weight) { compute_weight(num_chan, spec_data, mask_after_clipping_data, weight); set_weight_to_matrix(weight_matrix, irow, ipol, weight.at(var_index)); } } // end of polarization loop // output results of fitting if (num_apply_true == 0) continue; Array<Int> fpar_mtx(IPosition(2, num_pol, num_fpar_max), Array<Int>::uninitialized); set_matrix_for_bltable<size_t, Int>(num_pol, num_fpar_max, fpar_mtx_tmp, fpar_mtx); Array<Float> ffpar_mtx(IPosition(2, num_pol, num_ffpar_max), Array<Float>::uninitialized); set_matrix_for_bltable<double, Float>(num_pol, num_ffpar_max, ffpar_mtx_tmp, ffpar_mtx); Array<uInt> masklist_mtx(IPosition(2, num_pol, num_masklist_max), Array<uInt>::uninitialized); set_matrix_for_bltable<uInt, uInt>(num_pol, num_masklist_max, masklist_mtx_tmp, masklist_mtx); Array<Float> coeff_mtx(IPosition(2, num_pol, num_coeff_max), Array<Float>::uninitialized); set_matrix_for_bltable<double, Float>(num_pol, num_coeff_max, coeff_mtx_tmp, coeff_mtx); Matrix<uInt> masklist_mtx2 = masklist_mtx; Matrix<Bool> apply_mtx2 = apply_mtx; if (write_baseline_table) { bt->appenddata((uInt)scans[irow], (uInt)beams[irow], (uInt)antennas[irow], (uInt)data_spw[irow], 0, times[irow], apply_mtx, bltype_mtx, fpar_mtx, ffpar_mtx, masklist_mtx, coeff_mtx, rms_mtx, (uInt)num_chan, cthres_mtx, citer_mtx, uself_mtx, lfthres_mtx, lfavg_mtx, lfedge_mtx); } if (write_baseline_text) { for (size_t ipol = 0; ipol < num_pol; ++ipol) { if (apply_mtx2(ipol, 0) == false) continue; ofs_txt << "Scan" << '[' << (uInt)scans[irow] << ']' << ' ' << "Beam" << '[' << (uInt)beams[irow] << ']' << ' ' << "Spw" << '[' << (uInt)data_spw[irow] << ']' << ' ' << "Pol" << '[' << ipol << ']' << ' ' << "Time" << '[' << MVTime(times[irow]/ 24. / 3600.).string(MVTime::YMD, 8) << ']' << endl; ofs_txt << endl; ofs_txt << "Fitter range = " << '['; for (size_t imasklist = 0; imasklist < num_masklist_max/2; ++imasklist) { if (imasklist == 0) { ofs_txt << '[' << masklist_mtx2(ipol, 2 * imasklist) << ';' << masklist_mtx2(ipol, 2 * imasklist + 1) << ']'; } if (imasklist >= 1 && (0 != masklist_mtx2(ipol, 2 * imasklist) && 0 != masklist_mtx2(ipol, 2 * imasklist + 1))) { ofs_txt << ",[" << masklist_mtx2(ipol, 2 * imasklist) << ',' << masklist_mtx2(ipol, 2 * imasklist + 1) << ']'; } } ofs_txt << ']' << endl; ofs_txt << endl; Matrix<uInt> bltype_mtx2 = bltype_mtx[0][ipol]; Matrix<Int> fpar_mtx2 = fpar_mtx[0][ipol]; Matrix<Float> rms_mtx2 = rms_mtx[0][ipol]; string bltype_name; string blparam_name ="order"; if (bltype_mtx2(0, 0) == (uInt)0) { bltype_name = "poly"; } else if (bltype_mtx2(0, 0) == (uInt)1) { bltype_name = "chebyshev"; } else if (bltype_mtx2(0, 0) == (uInt)2) { blparam_name = "npiece"; bltype_name = "cspline"; } else if (bltype_mtx2(0, 0) == (uInt)3) { blparam_name = "nwave"; bltype_name = "sinusoid"; } ofs_txt << "Baseline parameters Function = " << bltype_name << " " << blparam_name << " = "; Matrix<Int> fpar_mtx3 = fpar_mtx; if (bltype_mtx2(0,0) == (uInt)3) { for (size_t num = 0; num < num_fpar_max; ++num) { ofs_txt << fpar_mtx3(ipol, num) << " "; } ofs_txt << endl; } else { ofs_txt << fpar_mtx2(0, 0) << endl; } ofs_txt << endl; ofs_txt << "Results of baseline fit" << endl; ofs_txt << endl; Matrix<Float> coeff_mtx2 = coeff_mtx; if (bltype_mtx2(0,0) == (uInt)0 || bltype_mtx2(0,0) == (uInt)1 || bltype_mtx2(0,0) == (uInt)2){ for (size_t icoeff = 0; icoeff < num_coeff_max; ++icoeff) { ofs_txt << "p" << icoeff << " = " << setprecision(8) << coeff_mtx2(ipol, icoeff) << " "; } } else if (bltype_mtx2(0,0) == (uInt)3) { size_t wn=0; string c_s ="s"; //if (blparam[0] == 0) { if (fpar_mtx3(ipol, wn) == 0) { ofs_txt << "c" << fpar_mtx3(ipol, wn) << " = " <<setw(13)<<left<< setprecision(8) << coeff_mtx2(ipol, 0) << " "; wn = 1; //for (size_t icoeff = 1; icoeff < num_coeff_max; ++icoeff) { for (size_t icoeff = 1; icoeff < coeff_mtx_tmp[ipol].size(); ++icoeff) { ofs_txt << c_s << fpar_mtx3(ipol, wn) << " = " <<setw(13)<<left<< setprecision(8) << coeff_mtx2(ipol, icoeff) << " "; c_s == "s" ? (c_s = "c") : (c_s = "s"); if (icoeff % 2 == 0) { ++wn; } } } else { wn = 0; //for (size_t icoeff = 0; icoeff < num_coeff_max; ++icoeff) { for (size_t icoeff = 0; icoeff < coeff_mtx_tmp[ipol].size(); ++icoeff) { ofs_txt << c_s << fpar_mtx3(ipol, wn) << " = " <<setw(13)<<left<< setprecision(8) << coeff_mtx2(ipol, icoeff) << " "; c_s == "s" ? (c_s = "c") : (c_s = "s"); if (icoeff % 2 != 0) { ++wn; } } } } ofs_txt << endl; ofs_txt << endl; ofs_txt << "rms = "; ofs_txt << setprecision(8) << rms_mtx2(0, 0) << endl; ofs_txt << endl; ofs_txt << "Number of clipped channels = " << final_mask_after_clipping[ipol] - final_mask[ipol] << endl; ofs_txt << endl; ofs_txt << "------------------------------------------------------" << endl; ofs_txt << endl; } } if (write_baseline_csv) { for (size_t ipol = 0; ipol < num_pol; ++ipol) { if (apply_mtx2(ipol, 0) == false) continue; ofs_csv << (uInt)scans[irow] << ',' << (uInt)beams[irow] << ',' << (uInt)data_spw[irow] << ',' << ipol << ',' << setprecision(12) << times[irow] << ','; ofs_csv << '['; for (size_t imasklist = 0; imasklist < num_masklist_max / 2; ++imasklist) { if (imasklist == 0) { ofs_csv << '[' << masklist_mtx2(ipol, 2 * imasklist) << ';' << masklist_mtx2(ipol, 2 * imasklist + 1) << ']'; } if (imasklist >= 1 && (0 != masklist_mtx2(ipol, 2 * imasklist) && 0 != masklist_mtx2(ipol, 2 * imasklist + 1))) { ofs_csv << ";[" << masklist_mtx2(ipol, 2 * imasklist) << ';' << masklist_mtx2(ipol, 2 * imasklist + 1) << ']'; } } ofs_csv << ']' << ','; Matrix<uInt> bltype_mtx2 = bltype_mtx[0][ipol]; string bltype_name; if (bltype_mtx2(0, 0) == (uInt)0) { bltype_name = "poly"; } else if (bltype_mtx2(0, 0) == (uInt)1) { bltype_name = "chebyshev"; } else if (bltype_mtx2(0, 0) == (uInt)2) { bltype_name = "cspline"; } else if (bltype_mtx2(0, 0) == (uInt)3) { bltype_name = "sinusoid"; } // TODO: revisit this line in CAS-13671 Matrix<Int> fpar_mtx2 = fpar_mtx; if (bltype_mtx2(0, 0) == (uInt)3) { ofs_csv << bltype_name.c_str() << ',' << fpar_mtx2(ipol, 0); for (size_t i = 1; i < num_fpar_max; ++i) { ofs_csv << ';' << fpar_mtx2(ipol, i); } ofs_csv << ','; } else { ofs_csv << bltype_name.c_str() << ',' << fpar_mtx2(ipol, 0) << ','; } Matrix<Float> coeff_mtx2 = coeff_mtx; if (bltype_mtx2(0, 0) == (uInt)3) { for (size_t icoeff = 0; icoeff < coeff_mtx_tmp[ipol].size(); ++icoeff) { ofs_csv << setprecision(8) << coeff_mtx2(ipol, icoeff) << ','; } } else { for (size_t icoeff = 0; icoeff < num_coeff_max; ++icoeff) { ofs_csv << setprecision(8) << coeff_mtx2(ipol, icoeff) << ','; } } Matrix<Float> rms_mtx2 = rms_mtx; ofs_csv << setprecision(8) << rms_mtx2(ipol, 0) << ','; ofs_csv << final_mask_after_clipping[ipol] - final_mask[ipol]; ofs_csv << endl; } } } // end of chunk row loop // write back data cube to VisBuffer if (do_subtract) { if (update_weight) { sdh_->fillCubeToOutputMs(vb, data_chunk, &flag_chunk, &weight_matrix); } else { sdh_->fillCubeToOutputMs(vb, data_chunk, &flag_chunk); } } } // end of vi loop } // end of chunk loop if (write_baseline_table) { bt->save(bloutputname_table); delete bt; } if (write_baseline_csv) { ofs_csv.close(); } if (write_baseline_text) { ofs_txt.close(); } finalize_process(); destroy_baseline_contexts(bl_contexts); //double tend = gettimeofday_sec(); //std::cout << "Elapsed time = " << (tend - tstart) << " sec." << std::endl; } //////////////////////////////////////////////////////////////////////// ///// Atcual processing functions //////////////////////////////////////////////////////////////////////// //Subtract baseline using normal or Chebyshev polynomials void SingleDishMS::subtractBaseline(string const& in_column_name, string const& out_ms_name, string const& out_bloutput_name, bool const& do_subtract, string const& in_spw, bool const& update_weight, string const& sigma_value, string const& blfunc, int const order, float const clip_threshold_sigma, int const num_fitting_max, bool const linefinding, float const threshold, int const avg_limit, int const minwidth, vector<int> const& edge) { vector<int> order_vect; order_vect.push_back(order); vector<int> blparam_exclude_dummy; blparam_exclude_dummy.clear(); LogIO os(_ORIGIN); os << "Fitting and subtracting polynomial baseline order = " << order << LogIO::POST; if (order < 0) { throw(AipsError("order must be positive or zero.")); } LIBSAKURA_SYMBOL(Status) status; LIBSAKURA_SYMBOL(LSQFitStatus) bl_status; std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> bl_contexts; bl_contexts.clear(); size_t bltype = BaselineType_kPolynomial; string blfunc_lower = blfunc; std::transform( blfunc_lower.begin(), blfunc_lower.end(), blfunc_lower.begin(), [](unsigned char c) {return std::tolower(c);} ); if (blfunc_lower == "chebyshev") { bltype = BaselineType_kChebyshev; } doSubtractBaseline(in_column_name, out_ms_name, out_bloutput_name, do_subtract, in_spw, update_weight, sigma_value, status, bl_contexts, bltype, order_vect, blparam_exclude_dummy, false, "", "", clip_threshold_sigma, num_fitting_max, linefinding, threshold, avg_limit, minwidth, edge, [&](LIBSAKURA_SYMBOL(LSQFitContextFloat) const *context, size_t const num_chan, std::vector<size_t> const &/*nwave*/, float *spec, bool const *mask, size_t const /*num_coeff*/, double *coeff, bool *mask_after_clipping, float *rms){ status = LIBSAKURA_SYMBOL(LSQFitPolynomialFloat)( context, static_cast<uint16_t>(order_vect[0]), num_chan, spec, mask, clip_threshold_sigma, num_fitting_max, order_vect[0] + 1, coeff, nullptr, nullptr, mask_after_clipping, rms, &bl_status); check_sakura_status("sakura_LSQFitPolynomialFloat", status); if (bl_status != LIBSAKURA_SYMBOL(LSQFitStatus_kOK)) { throw(AipsError("baseline fitting isn't successful.")); } }, [&](size_t ipol, std::vector<std::vector<double> > &ffpar_mtx_tmp, size_t &num_ffpar_max) { size_t num_ffpar = get_num_coeff_bloutput(bltype, 0, num_ffpar_max); ffpar_mtx_tmp[ipol].resize(num_ffpar); }, [&](LIBSAKURA_SYMBOL(LSQFitContextFloat) const *context, size_t const num_chan, std::vector<size_t> const &/*nwave*/, float *spec, size_t const /*num_coeff*/, double *coeff){ status = LIBSAKURA_SYMBOL(SubtractPolynomialFloat)( context, num_chan, spec, order_vect[0] + 1, coeff, spec); check_sakura_status("sakura_SubtractPolynomialFloat", status);}, [&](LIBSAKURA_SYMBOL(LSQFitContextFloat) const *context, size_t const num_chan, std::vector<size_t> const &/*nwave*/, size_t const /*num_coeff*/, float *spec, bool const *mask, bool *mask_after_clipping, float *rms){ status = LIBSAKURA_SYMBOL(LSQFitPolynomialFloat)( context, static_cast<uint16_t>(order_vect[0]), num_chan, spec, mask, clip_threshold_sigma, num_fitting_max, order_vect[0] + 1, nullptr, nullptr, spec, mask_after_clipping, rms, &bl_status); check_sakura_status("sakura_LSQFitPolynomialFloat", status); if (bl_status != LIBSAKURA_SYMBOL(LSQFitStatus_kOK)) { throw(AipsError("baseline fitting isn't successful.")); } }, os ); } //Subtract baseline using natural cubic spline void SingleDishMS::subtractBaselineCspline(string const& in_column_name, string const& out_ms_name, string const& out_bloutput_name, bool const& do_subtract, string const& in_spw, bool const& update_weight, string const& sigma_value, int const npiece, float const clip_threshold_sigma, int const num_fitting_max, bool const linefinding, float const threshold, int const avg_limit, int const minwidth, vector<int> const& edge) { vector<int> npiece_vect; npiece_vect.push_back(npiece); vector<int> blparam_exclude_dummy; blparam_exclude_dummy.clear(); LogIO os(_ORIGIN); os << "Fitting and subtracting cubic spline baseline npiece = " << npiece << LogIO::POST; if (npiece <= 0) { throw(AipsError("npiece must be positive.")); } LIBSAKURA_SYMBOL(Status) status; LIBSAKURA_SYMBOL(LSQFitStatus) bl_status; std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> bl_contexts; bl_contexts.clear(); size_t const bltype = BaselineType_kCubicSpline; SakuraAlignedArray<size_t> boundary(npiece+1); size_t *boundary_data = boundary.data(); doSubtractBaseline(in_column_name, out_ms_name, out_bloutput_name, do_subtract, in_spw, update_weight, sigma_value, status, bl_contexts, bltype, npiece_vect, blparam_exclude_dummy, false, "", "", clip_threshold_sigma, num_fitting_max, linefinding, threshold, avg_limit, minwidth, edge, [&](LIBSAKURA_SYMBOL(LSQFitContextFloat) const *context, size_t const num_chan, std::vector<size_t> const &/*nwave*/, float *spec, bool const *mask, size_t const /*num_coeff*/, double *coeff, bool *mask_after_clipping, float *rms) { status = LIBSAKURA_SYMBOL(LSQFitCubicSplineFloat)( context, static_cast<uint16_t>(npiece_vect[0]), num_chan, spec, mask, clip_threshold_sigma, num_fitting_max, reinterpret_cast<double (*)[4]>(coeff), nullptr, nullptr, mask_after_clipping, rms, boundary_data, &bl_status); check_sakura_status("sakura_LSQFitCubicSplineFloat", status); if (bl_status != LIBSAKURA_SYMBOL(LSQFitStatus_kOK)) { throw(AipsError("baseline fitting isn't successful.")); } }, [&](size_t ipol, std::vector<std::vector<double> > &ffpar_mtx_tmp, size_t &num_ffpar_max) { size_t num_ffpar = get_num_coeff_bloutput( bltype, npiece_vect[0], num_ffpar_max); ffpar_mtx_tmp[ipol].resize(num_ffpar); for (size_t ipiece = 0; ipiece < num_ffpar; ++ipiece) { ffpar_mtx_tmp[ipol][ipiece] = boundary_data[ipiece]; } }, [&](LIBSAKURA_SYMBOL(LSQFitContextFloat) const *context, size_t const num_chan, std::vector<size_t> const &/*nwave*/, float *spec, size_t const /*num_coeff*/, double *coeff) { status = LIBSAKURA_SYMBOL(SubtractCubicSplineFloat)( context, num_chan, spec, npiece_vect[0], reinterpret_cast<double (*)[4]>(coeff), boundary_data, spec); check_sakura_status("sakura_SubtractCubicSplineFloat", status);}, [&](LIBSAKURA_SYMBOL(LSQFitContextFloat) const *context, size_t const num_chan, std::vector<size_t> const &/*nwave*/, size_t const /*num_coeff*/, float *spec, bool const *mask, bool *mask_after_clipping, float *rms) { status = LIBSAKURA_SYMBOL(LSQFitCubicSplineFloat)( context, static_cast<uint16_t>(npiece_vect[0]), num_chan, spec, mask, clip_threshold_sigma, num_fitting_max, nullptr, nullptr, spec, mask_after_clipping, rms, boundary_data, &bl_status); check_sakura_status("sakura_LSQFitCubicSplineFloat", status); if (bl_status != LIBSAKURA_SYMBOL(LSQFitStatus_kOK)) { throw(AipsError("baseline fitting isn't successful.")); } }, os ); } void SingleDishMS::subtractBaselineSinusoid(string const& in_column_name, string const& out_ms_name, string const& out_bloutput_name, bool const& do_subtract, string const& in_spw, bool const& update_weight, string const& sigma_value, string const& addwn0, string const& rejwn0, bool const applyfft, string const fftmethod, string const fftthresh, float const clip_threshold_sigma, int const num_fitting_max, bool const linefinding, float const threshold, int const avg_limit, int const minwidth, vector<int> const& edge) { char const delim = ','; vector<int> addwn = string_to_list(addwn0, delim); vector<int> rejwn = string_to_list(rejwn0, delim); LogIO os(_ORIGIN); os << "Fitting and subtracting sinusoid baseline with wave numbers " << addwn0 << LogIO::POST; if (addwn.size() == 0) { throw(AipsError("addwn must contain at least one element.")); } LIBSAKURA_SYMBOL(Status) status; LIBSAKURA_SYMBOL(LSQFitStatus) bl_status; std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> bl_contexts; bl_contexts.clear(); LIBSAKURA_SYMBOL(LSQFitContextFloat) *context = nullptr; size_t bltype = BaselineType_kSinusoid; auto wn_ulimit_by_rejwn = [&](){ return ((rejwn.size() == 2) && (rejwn[1] == SinusoidWaveNumber_kUpperLimit)); }; auto par_spectrum_context = [&](){ return (applyfft && !wn_ulimit_by_rejwn()); }; auto prepare_context = [&](LIBSAKURA_SYMBOL(LSQFitContextFloat) const *context0, size_t const num_chan, std::vector<size_t> const &nwave){ if (par_spectrum_context()) { status = LIBSAKURA_SYMBOL(CreateLSQFitContextSinusoidFloat)( static_cast<uint16_t>(nwave[nwave.size()-1]), num_chan, &context); check_sakura_status("sakura_CreateLSQFitContextSinusoidFloat", status); } else { context = const_cast<LIBSAKURA_SYMBOL(LSQFitContextFloat) *>(context0); } }; auto clear_context = [&](){ if (par_spectrum_context()) { status = LIBSAKURA_SYMBOL(DestroyLSQFitContextFloat)(context); check_sakura_status("sakura_DestoyBaselineContextFloat", status); context = nullptr; } }; doSubtractBaseline(in_column_name, out_ms_name, out_bloutput_name, do_subtract, in_spw, update_weight, sigma_value, status, bl_contexts, bltype, addwn, rejwn, applyfft, fftmethod, fftthresh, clip_threshold_sigma, num_fitting_max, linefinding, threshold, avg_limit, minwidth, edge, [&](LIBSAKURA_SYMBOL(LSQFitContextFloat) const *context0, size_t const num_chan, std::vector<size_t> const &nwave, float *spec, bool const *mask, size_t const num_coeff, double *coeff, bool *mask_after_clipping, float *rms) { prepare_context(context0, num_chan, nwave); status = LIBSAKURA_SYMBOL(LSQFitSinusoidFloat)( context, nwave.size(), &nwave[0], num_chan, spec, mask, clip_threshold_sigma, num_fitting_max, num_coeff, coeff, nullptr, nullptr, mask_after_clipping, rms, &bl_status); check_sakura_status("sakura_LSQFitSinusoidFloat", status); check_baseline_status(bl_status); }, [&](size_t ipol, std::vector<std::vector<double> > &ffpar_mtx_tmp, size_t &num_ffpar_max) { size_t num_ffpar = get_num_coeff_bloutput(bltype, addwn.size(), num_ffpar_max); ffpar_mtx_tmp[ipol].resize(num_ffpar); }, [&](LIBSAKURA_SYMBOL(LSQFitContextFloat) const *context0, size_t const num_chan, std::vector<size_t> const &nwave, float *spec, size_t num_coeff, double *coeff) { if (!par_spectrum_context()) { context = const_cast<LIBSAKURA_SYMBOL(LSQFitContextFloat) *>(context0); } status = LIBSAKURA_SYMBOL(SubtractSinusoidFloat)( context, num_chan, spec, nwave.size(), &nwave[0], num_coeff, coeff, spec); check_sakura_status("sakura_SubtractSinusoidFloat", status); clear_context(); }, [&](LIBSAKURA_SYMBOL(LSQFitContextFloat) const *context0, size_t const num_chan, std::vector<size_t> const &nwave, size_t const num_coeff, float *spec, bool const *mask, bool *mask_after_clipping, float *rms) { prepare_context(context0, num_chan, nwave); status = LIBSAKURA_SYMBOL(LSQFitSinusoidFloat)( context, nwave.size(), &nwave[0], num_chan, spec, mask, clip_threshold_sigma, num_fitting_max, num_coeff, nullptr, nullptr, spec, mask_after_clipping, rms, &bl_status); check_sakura_status("sakura_LSQFitSinusoidFloat", status); check_baseline_status(bl_status); clear_context(); }, os ); } // Apply baseline table to MS void SingleDishMS::applyBaselineTable(string const& in_column_name, string const& in_bltable_name, string const& in_spw, bool const& update_weight, string const& sigma_value, string const& out_ms_name) { LogIO os(_ORIGIN); os << "Apply baseline table " << in_bltable_name << " to MS. " << LogIO::POST; if (in_bltable_name == "") { throw(AipsError("baseline table is not given.")); } // parse fitting parameters in the text file BLTableParser parser(in_bltable_name); std::vector<size_t> baseline_types = parser.get_function_types(); map<size_t const, uint16_t> max_orders; for (size_t i = 0; i < baseline_types.size(); ++i) { max_orders[baseline_types[i]] = parser.get_max_order(baseline_types[i]); } { //DEBUG output os << LogIO::DEBUG1 << "spw ID = " << in_spw << LogIO::POST; os << "Baseline Types = " << baseline_types << LogIO::POST; os << "Max Orders:" << LogIO::POST; map<size_t const, uint16_t>::iterator iter = max_orders.begin(); while (iter != max_orders.end()) { os << LogIO::DEBUG1 << "- type " << (*iter).first << ": " << (*iter).second << LogIO::POST; ++iter; } } Block<Int> columns(1); columns[0] = MS::DATA_DESC_ID; // (CAS-9918, 2017/4/27 WK) //columns[0] = MS::TIME; prepare_for_process(in_column_name, out_ms_name, columns, false); vi::VisibilityIterator2 *vi = sdh_->getVisIter(); vi->setRowBlocking(kNRowBlocking); vi::VisBuffer2 *vb = vi->getVisBuffer(); Vector<Int> recspw; Matrix<Int> recchan; Vector<size_t> nchan; Vector<Vector<Bool> > in_mask; Vector<bool> nchan_set; parse_spw(in_spw, recspw, recchan, nchan, in_mask, nchan_set); // Baseline Contexts reservoir // key: BaselineType // value: a vector of Sakura_BaselineContextFloat for various nchans Vector<size_t> ctx_indices(nchan.nelements(), 0ul); map<size_t const, std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> > context_reservoir; map<size_t const, uint16_t>::iterator iter = max_orders.begin(); while (iter != max_orders.end()) { context_reservoir[(*iter).first] = std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *>(); ++iter; } LIBSAKURA_SYMBOL(Status) status; std::vector<float> weight(WeightIndex_kNum); size_t const var_index = (sigma_value == "stddev") ? WeightIndex_kStddev : WeightIndex_kRms; for (vi->originChunks(); vi->moreChunks(); vi->nextChunk()) { for (vi->origin(); vi->more(); vi->next()) { Vector<Int> scans = vb->scan(); Vector<Double> times = vb->time(); Vector<Double> intervals = vb->timeInterval(); Vector<Int> beams = vb->feed1(); Vector<Int> antennas = vb->antenna1(); Vector<Int> data_spw = vb->spectralWindows(); size_t const num_chan = static_cast<size_t>(vb->nChannels()); size_t const num_pol = static_cast<size_t>(vb->nCorrelations()); size_t const num_row = static_cast<size_t>(vb->nRows()); Cube<Float> data_chunk(num_pol, num_chan, num_row); SakuraAlignedArray<float> spec(num_chan); float *spec_data = spec.data(); Cube<Bool> flag_chunk(num_pol, num_chan, num_row); SakuraAlignedArray<bool> flag(num_chan); bool *flag_data = flag.data(); SakuraAlignedArray<bool> mask(num_chan); bool *mask_data = mask.data(); Matrix<Float> weight_matrix(num_pol, num_row, Array<Float>::uninitialized); bool new_nchan = false; get_nchan_and_mask(recspw, data_spw, recchan, num_chan, nchan, in_mask, nchan_set, new_nchan); if (new_nchan) { map<size_t const, uint16_t>::iterator iter = max_orders.begin(); while (iter != max_orders.end()) { get_baseline_context((*iter).first, (*iter).second, num_chan, nchan, nchan_set, ctx_indices, context_reservoir[(*iter).first]); ++iter; } } // get data/flag cubes (npol*nchan*nrow) from VisBuffer get_data_cube_float(*vb, data_chunk); get_flag_cube(*vb, flag_chunk); // get weight matrix (npol*nrow) from VisBuffer if (update_weight) { get_weight_matrix(*vb, weight_matrix); } // loop over MS rows for (size_t irow = 0; irow < num_row; ++irow) { size_t idx = 0; for (size_t ispw = 0; ispw < recspw.nelements(); ++ispw) { if (data_spw[irow] == recspw[ispw]) { idx = ispw; break; } } // find a baseline table row (index) corresponding to this MS row size_t idx_fit_param; if (!parser.GetFitParameterIdx(times[irow], intervals[irow], scans[irow], beams[irow], antennas[irow], data_spw[irow], idx_fit_param)) { for (size_t ipol = 0; ipol < num_pol; ++ipol) { flag_spectrum_in_cube(flag_chunk, irow, ipol); //flag } continue; } // loop over polarization for (size_t ipol = 0; ipol < num_pol; ++ipol) { bool apply; Vector<double> coeff; Vector<size_t> boundary; std::vector<bool> mask_bltable; BLParameterSet fit_param; parser.GetFitParameterByIdx(idx_fit_param, ipol, apply, coeff, boundary, mask_bltable, fit_param); if (!apply) { flag_spectrum_in_cube(flag_chunk, irow, ipol); //flag continue; } // get a channel mask from data cube // (note that the variable 'mask' is flag in the next line // actually, then it will be converted to real mask when // taking AND with user-given mask info. this is just for // saving memory usage...) get_flag_from_cube(flag_chunk, irow, ipol, num_chan, flag_data); // skip spectrum if all channels flagged if (allchannels_flagged(num_chan, flag_data)) { continue; } // get a spectrum from data cube get_spectrum_from_cube(data_chunk, irow, ipol, num_chan, spec_data); // actual execution of single spectrum map<size_t const, std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> >::iterator iter = context_reservoir.find(fit_param.baseline_type); if (iter == context_reservoir.end()) throw(AipsError("Invalid baseline type detected!")); LIBSAKURA_SYMBOL(LSQFitContextFloat) * context = (*iter).second[ctx_indices[idx]]; //cout << "Got context for type " << (*iter).first << ": idx=" << ctx_indices[idx] << endl; SakuraAlignedArray<double> coeff_storage(coeff); double *coeff_data = coeff_storage.data(); SakuraAlignedArray<size_t> boundary_storage(boundary); size_t *boundary_data = boundary_storage.data(); string subtract_funcname; switch (static_cast<size_t>(fit_param.baseline_type)) { case BaselineType_kPolynomial: case BaselineType_kChebyshev: status = LIBSAKURA_SYMBOL(SubtractPolynomialFloat)( context, num_chan, spec_data, coeff.size(), coeff_data, spec_data); subtract_funcname = "sakura_SubtractPolynomialFloat"; break; case BaselineType_kCubicSpline: status = LIBSAKURA_SYMBOL(SubtractCubicSplineFloat)( context, num_chan, spec_data, boundary.size()-1, reinterpret_cast<double (*)[4]>(coeff_data), boundary_data, spec_data); subtract_funcname = "sakura_SubtractCubicSplineFloat"; break; case BaselineType_kSinusoid: status = LIBSAKURA_SYMBOL(SubtractSinusoidFloat)( context, num_chan, spec_data, fit_param.nwave.size(), &fit_param.nwave[0], coeff.size(), coeff_data, spec_data); subtract_funcname = "sakura_SubtractSinusoidFloat"; break; default: throw(AipsError("Unsupported baseline type.")); } check_sakura_status(subtract_funcname, status); // set back a spectrum to data cube set_spectrum_to_cube(data_chunk, irow, ipol, num_chan, spec_data); if (update_weight) { // convert flag to mask by taking logical NOT of flag // and then operate logical AND with in_mask and with mask from bltable for (size_t ichan = 0; ichan < num_chan; ++ichan) { mask_data[ichan] = in_mask[idx][ichan] && (!(flag_data[ichan])) && mask_bltable[ichan]; } compute_weight(num_chan, spec_data, mask_data, weight); set_weight_to_matrix(weight_matrix, irow, ipol, weight.at(var_index)); } } // end of polarization loop } // end of chunk row loop // write back data and flag cube to VisBuffer if (update_weight) { sdh_->fillCubeToOutputMs(vb, data_chunk, &flag_chunk, &weight_matrix); } else { sdh_->fillCubeToOutputMs(vb, data_chunk, &flag_chunk); } } // end of vi loop } // end of chunk loop finalize_process(); // destroy baseline contexts map<size_t const, std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> >::iterator ctxiter = context_reservoir.begin(); while (ctxiter != context_reservoir.end()) { destroy_baseline_contexts (context_reservoir[(*ctxiter).first]); ++ctxiter; } } // Fit line profile void SingleDishMS::fitLine(string const& in_column_name, string const& in_spw, string const& /* in_pol */, string const& fitfunc, string const& in_nfit, bool const linefinding, float const threshold, int const avg_limit, int const minwidth, vector<int> const& edge, string const& tempfile_name, string const& temp_out_ms_name) { // in_column = [FLOAT_DATA|DATA|CORRECTED_DATA] // no iteration is necessary for the processing. // procedure // 1. iterate over MS // 2. pick single spectrum from in_column (this is not actually necessary for simple scaling but for exibision purpose) // 3. fit Gaussian or Lorentzian profile to each spectrum // 4. write fitting results to outfile LogIO os(_ORIGIN); os << "Fitting line profile with " << fitfunc << LogIO::POST; if (file_exists(temp_out_ms_name)) { throw(AipsError("temporary ms file unexpectedly exists.")); } if (file_exists(tempfile_name)) { throw(AipsError("temporary file unexpectedly exists.")); } Block<Int> columns(1); columns[0] = MS::DATA_DESC_ID; prepare_for_process(in_column_name, temp_out_ms_name, columns, false); vi::VisibilityIterator2 *vi = sdh_->getVisIter(); vi->setRowBlocking(kNRowBlocking); vi::VisBuffer2 *vb = vi->getVisBuffer(); ofstream ofs(tempfile_name); Vector<Int> recspw; Matrix<Int> recchan; Vector<size_t> nchan; Vector<Vector<Bool> > in_mask; Vector<bool> nchan_set; parse_spw(in_spw, recspw, recchan, nchan, in_mask, nchan_set); std::vector<size_t> nfit; if (linefinding) { os << "Defining line ranges using line finder. nfit will be ignored." << LogIO::POST; } else { std::vector<string> nfit_s = split_string(in_nfit, ','); nfit.resize(nfit_s.size()); for (size_t i = 0; i < nfit_s.size(); ++i) { nfit[i] = std::stoi(nfit_s[i]); } } size_t num_spec = 0; size_t num_noline = 0; for (vi->originChunks(); vi->moreChunks(); vi->nextChunk()) { for (vi->origin(); vi->more(); vi->next()) { Vector<Int> scans = vb->scan(); Vector<Double> times = vb->time(); Vector<Int> beams = vb->feed1(); Vector<Int> antennas = vb->antenna1(); Vector<Int> data_spw = vb->spectralWindows(); size_t const num_chan = static_cast<size_t>(vb->nChannels()); size_t const num_pol = static_cast<size_t>(vb->nCorrelations()); size_t const num_row = static_cast<size_t>(vb->nRows()); Cube<Float> data_chunk(num_pol, num_chan, num_row); SakuraAlignedArray<float> spec(num_chan); Cube<Bool> flag_chunk(num_pol, num_chan, num_row); SakuraAlignedArray<bool> mask(num_chan); // CAUTION!!! // data() method must be used with special care!!! float *spec_data = spec.data(); bool *mask_data = mask.data(); bool new_nchan = false; get_nchan_and_mask(recspw, data_spw, recchan, num_chan, nchan, in_mask, nchan_set, new_nchan); // get data/flag cubes (npol*nchan*nrow) from VisBuffer get_data_cube_float(*vb, data_chunk); get_flag_cube(*vb, flag_chunk); // loop over MS rows for (size_t irow = 0; irow < num_row; ++irow) { size_t idx = 0; for (size_t ispw = 0; ispw < recspw.nelements(); ++ispw) { if (data_spw[irow] == recspw[ispw]) { idx = ispw; break; } } std::vector<size_t> fitrange_start; fitrange_start.clear(); std::vector<size_t> fitrange_end; fitrange_end.clear(); for (size_t i = 0; i < recchan.nrow(); ++i) { if (recchan.row(i)(0) == data_spw[irow]) { fitrange_start.push_back(recchan.row(i)(1)); fitrange_end.push_back(recchan.row(i)(2)); } } if (!linefinding && nfit.size() != fitrange_start.size()) { throw(AipsError( "the number of elements of nfit and fitranges specified in spw must be identical.")); } // loop over polarization for (size_t ipol = 0; ipol < num_pol; ++ipol) { // get a channel mask from data cube // (note that the variable 'mask' is flag in the next line // actually, then it will be converted to real mask when // taking AND with user-given mask info. this is just for // saving memory usage...) get_flag_from_cube(flag_chunk, irow, ipol, num_chan, mask_data); // skip spectrum if all channels flagged if (allchannels_flagged(num_chan, mask_data)) { continue; } ++num_spec; // convert flag to mask by taking logical NOT of flag // and then operate logical AND with in_mask for (size_t ichan = 0; ichan < num_chan; ++ichan) { mask_data[ichan] = in_mask[idx][ichan] && (!(mask_data[ichan])); } // get a spectrum from data cube get_spectrum_from_cube(data_chunk, irow, ipol, num_chan, spec_data); // line finding. get fit mask (invert=false) if (linefinding) { list<pair<size_t, size_t>> line_ranges = findLineAndGetRanges(num_chan, spec_data, mask_data, threshold, avg_limit, minwidth, edge, false); if (line_ranges.size()==0) { ++num_noline; continue; } size_t nline = line_ranges.size(); fitrange_start.resize(nline); fitrange_end.resize(nline); nfit.resize(nline); auto range=line_ranges.begin(); for (size_t iline=0; iline<nline; ++iline){ fitrange_start[iline] = (*range).first; fitrange_end[iline] = (*range).second; nfit[iline] = 1; ++range; } } Vector<Float> x_; x_.resize(num_chan); Vector<Float> y_; y_.resize(num_chan); Vector<Bool> m_; m_.resize(num_chan); for (size_t ichan = 0; ichan < num_chan; ++ichan) { x_[ichan] = static_cast<Float>(ichan); y_[ichan] = spec_data[ichan]; } Vector<Float> parameters_; Vector<Float> error_; PtrBlock<Function<Float>*> funcs_; std::vector<std::string> funcnames_; std::vector<int> funccomponents_; std::string expr; if (fitfunc == "gaussian") { expr = "gauss"; } else if (fitfunc == "lorentzian") { expr = "lorentz"; } bool any_converged = false; for (size_t ifit = 0; ifit < nfit.size(); ++ifit) { if (nfit[ifit] == 0) continue; if (0 < ifit) ofs << ":"; //extract spec/mask within fitrange for (size_t ichan = 0; ichan < num_chan; ++ichan) { if ((fitrange_start[ifit] <= ichan) && (ichan <= fitrange_end[ifit])) { m_[ichan] = mask_data[ichan]; } else { m_[ichan] = false; } } //initial guesss Vector<Float> peak; Vector<Float> cent; Vector<Float> fwhm; peak.resize(nfit[ifit]); cent.resize(nfit[ifit]); fwhm.resize(nfit[ifit]); if (nfit[ifit] == 1) { Float sum = 0.0; Float max_spec = fabs(y_[fitrange_start[ifit]]); Float max_spec_x = x_[fitrange_start[ifit]]; bool is_positive = true; for (size_t ichan = fitrange_start[ifit]; ichan <= fitrange_end[ifit]; ++ichan) { sum += y_[ichan]; if (max_spec < fabs(y_[ichan])) { max_spec = fabs(y_[ichan]); max_spec_x = x_[ichan]; is_positive = (fabs(y_[ichan]) == y_[ichan]); } } peak[0] = max_spec * (is_positive ? 1 : -1); cent[0] = max_spec_x; fwhm[0] = fabs(sum / max_spec * 0.7); } else { size_t x_start = fitrange_start[ifit]; size_t x_width = (fitrange_end[ifit] - fitrange_start[ifit]) / nfit[ifit]; size_t x_end = x_start + x_width; for (size_t icomp = 0; icomp < nfit[ifit]; ++icomp) { if (icomp == nfit[ifit] - 1) { x_end = fitrange_end[ifit] + 1; } Float sum = 0.0; Float max_spec = fabs(y_[x_start]); Float max_spec_x = x_[x_start]; bool is_positive = true; for (size_t ichan = x_start; ichan < x_end; ++ichan) { sum += y_[ichan]; if (max_spec < fabs(y_[ichan])) { max_spec = fabs(y_[ichan]); max_spec_x = x_[ichan]; is_positive = (fabs(y_[ichan]) == y_[ichan]); } } peak[icomp] = max_spec * (is_positive ? 1 : -1); cent[icomp] = max_spec_x; fwhm[icomp] = fabs(sum / max_spec * 0.7); x_start += x_width; x_end += x_width; } } //fitter setup funcs_.resize(nfit[ifit]); funcnames_.clear(); funccomponents_.clear(); for (size_t icomp = 0; icomp < funcs_.nelements(); ++icomp) { if (expr == "gauss") { funcs_[icomp] = new Gaussian1D<Float>(); } else if (expr == "lorentz") { funcs_[icomp] = new Lorentzian1D<Float>(); } (funcs_[icomp]->parameters())[0] = peak[icomp]; //initial guess (peak) (funcs_[icomp]->parameters())[1] = cent[icomp]; //initial guess (centre) (funcs_[icomp]->parameters())[2] = fwhm[icomp]; //initial guess (fwhm) funcnames_.push_back(expr); funccomponents_.push_back(3); } //actual fitting NonLinearFitLM<Float> fitter; CompoundFunction<Float> func; for (size_t icomp = 0; icomp < funcs_.nelements(); ++icomp) { func.addFunction(*funcs_[icomp]); } fitter.setFunction(func); fitter.setMaxIter(50 + 10 * funcs_.nelements()); fitter.setCriteria(0.001); // Convergence criterium parameters_.resize(); parameters_ = fitter.fit(x_, y_, &m_); any_converged |= fitter.converged(); // if (!fitter.converged()) { // throw(AipsError("Failed in fitting. Fitter did not converge.")); // } error_.resize(); error_ = fitter.errors(); //write best-fit parameters to tempfile/outfile for (size_t icomp = 0; icomp < funcs_.nelements(); ++icomp) { if (0 < icomp) ofs << ":"; size_t offset = 3 * icomp; ofs.precision(4); ofs.setf(ios::fixed); ofs << scans[irow] << "," // scanID << times[irow] << "," // time << antennas[irow] << "," // antennaID << beams[irow] << "," // beamID << data_spw[irow] << "," // spwID << ipol << ","; // polID ofs.precision(8); ofs << parameters_[offset + 1] << "," << error_[offset + 1] << "," // cent << parameters_[offset + 0] << "," << error_[offset + 0] << "," // peak << parameters_[offset + 2] << "," << error_[offset + 2]; // fwhm } } //end of nfit loop ofs << "\n"; // count up spectra w/o any line fit if (!any_converged) ++num_noline; } //end of polarization loop } // end of MS row loop } //end of vi loop } //end of chunk loop if (num_noline==num_spec) { os << LogIO::WARN << "Fitter did not converge on any fit components." << LogIO::POST; } else if (num_noline > 0) { os << "No convergence for fitting to " << num_noline << " out of " << num_spec << " spectra" << LogIO::POST; } finalize_process(); ofs.close(); } //Subtract baseline by per spectrum fitting parameters void SingleDishMS::subtractBaselineVariable(string const& in_column_name, string const& out_ms_name, string const& out_bloutput_name, bool const& do_subtract, string const& in_spw, bool const& update_weight, string const& sigma_value, string const& param_file, bool const& verbose) { LogIO os(_ORIGIN); os << "Fitting and subtracting baseline using parameters in file " << param_file << LogIO::POST; Block<Int> columns(1); columns[0] = MS::DATA_DESC_ID; prepare_for_process(in_column_name, out_ms_name, columns, false); vi::VisibilityIterator2 *vi = sdh_->getVisIter(); vi->setRowBlocking(kNRowBlocking); vi::VisBuffer2 *vb = vi->getVisBuffer(); split_bloutputname(out_bloutput_name); bool write_baseline_csv = (bloutputname_csv != ""); bool write_baseline_text = (bloutputname_text != ""); bool write_baseline_table = (bloutputname_table != ""); ofstream ofs_csv; ofstream ofs_txt; BaselineTable *bt = 0; if (write_baseline_csv) { ofs_csv.open(bloutputname_csv.c_str()); } if (write_baseline_text) { ofs_txt.open(bloutputname_text.c_str(), std::ios::app); } if (write_baseline_table) { bt = new BaselineTable(vi->ms()); } Vector<Int> recspw; Matrix<Int> recchan; Vector<size_t> nchan; Vector<Vector<Bool> > in_mask; Vector<bool> nchan_set; parse_spw(in_spw, recspw, recchan, nchan, in_mask, nchan_set); // parse fitting parameters in the text file BLParameterParser parser(param_file); std::vector<size_t> baseline_types = parser.get_function_types(); /* max_orders: { baseline type as from enum, or poly/chebyshev: order or cspline: npiece or sinusoid: nwave.size() } Note: the biggest one of each? */ map<size_t const, uint16_t> max_orders; for (size_t i = 0; i < baseline_types.size(); ++i) { max_orders[baseline_types[i]] = parser.get_max_order(baseline_types[i]); } { //DEBUG ouput os << "Baseline Types = " << baseline_types << LogIO::POST; os << "Max Orders:" << LogIO::POST; map<size_t const, uint16_t>::iterator iter = max_orders.begin(); while (iter != max_orders.end()) { os << LogIO::DEBUG1 << "- type " << (*iter).first << ": " << (*iter).second << LogIO::POST; ++iter; } } // Baseline Contexts reservoir // key: Sakura_BaselineType enum, // value: a vector of Sakura_BaselineContextFloat for various nchans map<size_t const, std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> > context_reservoir; { map<size_t const, uint16_t>::iterator iter = max_orders.begin(); while (iter != max_orders.end()) { context_reservoir[(*iter).first] = std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *>(); ++iter; } } Vector<size_t> ctx_indices(recspw.size(), 0ul); //stores the number of channels of corresponding elements in contexts list. // WORKAROUND for absense of the way to get num_bases_data in context. vector<size_t> ctx_nchans; LIBSAKURA_SYMBOL(Status) status; LIBSAKURA_SYMBOL(LSQFitStatus) bl_status; std::vector<float> weight(WeightIndex_kNum); size_t const var_index = (sigma_value == "stddev") ? WeightIndex_kStddev : WeightIndex_kRms; for (vi->originChunks(); vi->moreChunks(); vi->nextChunk()) { for (vi->origin(); vi->more(); vi->next()) { Vector<Int> scans = vb->scan(); Vector<Double> times = vb->time(); Vector<Int> beams = vb->feed1(); Vector<Int> antennas = vb->antenna1(); Vector<Int> data_spw = vb->spectralWindows(); size_t const num_chan = static_cast<size_t>(vb->nChannels()); size_t const num_pol = static_cast<size_t>(vb->nCorrelations()); size_t const num_row = static_cast<size_t>(vb->nRows()); auto orig_rows = vb->rowIds(); Cube<Float> data_chunk(num_pol, num_chan, num_row); SakuraAlignedArray<float> spec(num_chan); Cube<Bool> flag_chunk(num_pol, num_chan, num_row); SakuraAlignedArray<bool> flag(num_chan); SakuraAlignedArray<bool> mask(num_chan); SakuraAlignedArray<bool> mask_after_clipping(num_chan); // CAUTION!!! // data() method must be used with special care!!! float *spec_data = spec.data(); bool *flag_data = flag.data(); bool *mask_data = mask.data(); bool *mask_after_clipping_data = mask_after_clipping.data(); Matrix<Float> weight_matrix(num_pol, num_row, Array<Float>::uninitialized); uInt final_mask[num_pol]; uInt final_mask_after_clipping[num_pol]; final_mask[0] = 0; final_mask[1] = 0; final_mask_after_clipping[0] = 0; final_mask_after_clipping[1] = 0; bool new_nchan = false; get_nchan_and_mask(recspw, data_spw, recchan, num_chan, nchan, in_mask, nchan_set, new_nchan); // check if context should be created once per chunk // in the first actual excution of baseline. bool check_context = true; // get data/flag cubes (npol*nchan*nrow) from VisBuffer get_data_cube_float(*vb, data_chunk); get_flag_cube(*vb, flag_chunk); // get weight matrix (npol*nrow) from VisBuffer if (update_weight) { get_weight_matrix(*vb, weight_matrix); } // loop over MS rows for (size_t irow = 0; irow < num_row; ++irow) { size_t idx = 0; for (size_t ispw = 0; ispw < recspw.nelements(); ++ispw) { if (data_spw[irow] == recspw[ispw]) { idx = ispw; break; } } //prepare varables for writing baseline table Array<Bool> apply_mtx(IPosition(2, num_pol, 1), true); Array<uInt> bltype_mtx(IPosition(2, num_pol, 1)); Array<Int> fpar_mtx(IPosition(2, num_pol, 1)); std::vector<std::vector<double> > ffpar_mtx_tmp(num_pol); std::vector<std::vector<uInt> > masklist_mtx_tmp(num_pol); std::vector<std::vector<double> > coeff_mtx_tmp(num_pol); Array<Float> rms_mtx(IPosition(2, num_pol, 1)); Array<Float> cthres_mtx(IPosition(2, num_pol, 1)); Array<uInt> citer_mtx(IPosition(2, num_pol, 1)); Array<Bool> uself_mtx(IPosition(2, num_pol, 1)); Array<Float> lfthres_mtx(IPosition(2, num_pol, 1)); Array<uInt> lfavg_mtx(IPosition(2, num_pol, 1)); Array<uInt> lfedge_mtx(IPosition(2, num_pol, 2)); size_t num_apply_true = 0; size_t num_ffpar_max = 0; size_t num_masklist_max = 0; size_t num_coeff_max = 0; std::vector<size_t> numcoeff(num_pol); // loop over polarization for (size_t ipol = 0; ipol < num_pol; ++ipol) { // get a channel mask from data cube // (note that the variable 'mask' is flag in the next line // actually, then it will be converted to real mask when // taking AND with user-given mask info. this is just for // saving memory usage...) get_flag_from_cube(flag_chunk, irow, ipol, num_chan, flag_data); // skip spectrum if all channels flagged if (allchannels_flagged(num_chan, flag_data)) { os << LogIO::DEBUG1 << "Row " << orig_rows[irow] << ", Pol " << ipol << ": All channels flagged. Skipping." << LogIO::POST; apply_mtx[0][ipol] = false; continue; } // convert flag to mask by taking logical NOT of flag // and then operate logical AND with in_mask for (size_t ichan = 0; ichan < num_chan; ++ichan) { mask_data[ichan] = in_mask[idx][ichan] && (!(flag_data[ichan])); } // get fitting parameter BLParameterSet fit_param; if (!parser.GetFitParameter(orig_rows[irow], ipol, fit_param)) { //no fit requrested flag_spectrum_in_cube(flag_chunk, irow, ipol); os << LogIO::DEBUG1 << "Row " << orig_rows[irow] << ", Pol " << ipol << ": Fit not requested. Skipping." << LogIO::POST; apply_mtx[0][ipol] = false; continue; } if (verbose) { os << "Fitting Parameter" << LogIO::POST; os << "[ROW " << orig_rows[irow] << " (nchan " << num_chan << ")" << ", POL" << ipol << "]" << LogIO::POST; fit_param.PrintSummary(); } // Create contexts when actually subtract baseine for the first time (if not yet exist) if (check_context) { // Generate context for all necessary baseline types map<size_t const, uint16_t>::iterator iter = max_orders.begin(); while (iter != max_orders.end()) { get_baseline_context((*iter).first, (*iter).second, num_chan, idx, ctx_indices, ctx_nchans, context_reservoir[(*iter).first]); ++iter; } check_context = false; } // get mask from BLParameterset and create composit mask if (fit_param.baseline_mask != "") { stringstream local_spw; local_spw << data_spw[irow] << ":" << fit_param.baseline_mask; Record selrec = sdh_->getSelRec(local_spw.str()); Matrix<Int> local_rec_chan = selrec.asArrayInt("channel"); Vector<Bool> local_mask(num_chan, false); get_mask_from_rec(data_spw[irow], local_rec_chan, local_mask, false); for (size_t ichan = 0; ichan < num_chan; ++ichan) { mask_data[ichan] = mask_data[ichan] && local_mask[ichan]; } } // check for composit mask and flag if no valid channel to fit if (NValidMask(num_chan, mask_data) == 0) { flag_spectrum_in_cube(flag_chunk, irow, ipol); apply_mtx[0][ipol] = false; os << LogIO::DEBUG1 << "Row " << orig_rows[irow] << ", Pol " << ipol << ": No valid channel to fit. Skipping" << LogIO::POST; continue; } // get a spectrum from data cube get_spectrum_from_cube(data_chunk, irow, ipol, num_chan, spec_data); // get baseline context map<size_t const, std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> >::iterator iter = context_reservoir.find(fit_param.baseline_type); if (iter == context_reservoir.end()) { throw(AipsError("Invalid baseline type detected!")); } LIBSAKURA_SYMBOL(LSQFitContextFloat) * context = (*iter).second[ctx_indices[idx]]; // Number of coefficients to fit this spectrum size_t num_coeff; size_t bltype = static_cast<size_t>(fit_param.baseline_type); switch (bltype) { case BaselineType_kPolynomial: case BaselineType_kChebyshev: status = LIBSAKURA_SYMBOL(GetNumberOfCoefficientsFloat)(context, fit_param.order, &num_coeff); check_sakura_status("sakura_GetNumberOfCoefficientsFloat", status); break; case BaselineType_kCubicSpline: num_coeff = 4 * fit_param.npiece; break; case BaselineType_kSinusoid: /* From sakuralib docs: The number of elements in the array coeff. If coeff is not null pointer, it must be ( num_nwave*2-1 ) or ( num_nwave*2 ) in cases nwave contains zero or not, respectively, and must not exceed num_data, while the value is not checked when coeff is null pointer. */ if (fit_param.nwave[0] == 0) num_coeff = fit_param.nwave.size() * 2 - 1; else num_coeff = fit_param.nwave.size() * 2; break; default: throw(AipsError("Unsupported baseline type.")); } numcoeff[ipol] = num_coeff; // Final check of the valid number of channels size_t num_min = (bltype == BaselineType_kCubicSpline) ? fit_param.npiece + 3 : num_coeff; if (NValidMask(num_chan, mask_data) < num_min) { flag_spectrum_in_cube(flag_chunk, irow, ipol); apply_mtx[0][ipol] = false; os << LogIO::WARN << "Too few valid channels to fit. Skipping Antenna " << antennas[irow] << ", Beam " << beams[irow] << ", SPW " << data_spw[irow] << ", Pol " << ipol << ", Time " << MVTime(times[irow] / 24. / 3600.).string(MVTime::YMD, 8) << LogIO::POST; continue; } // actual execution of single spectrum float rms; size_t num_boundary = 0; if (bltype == BaselineType_kCubicSpline) { num_boundary = fit_param.npiece+1; } SakuraAlignedArray<size_t> boundary(num_boundary); size_t *boundary_data = boundary.data(); if (write_baseline_text || write_baseline_csv || write_baseline_table) { num_apply_true++; bltype_mtx[0][ipol] = (uInt)fit_param.baseline_type; Vector<Int> fpar_tmp(1,0); switch (bltype) { case BaselineType_kPolynomial: case BaselineType_kChebyshev: fpar_tmp = (Int)fit_param.order; break; case BaselineType_kCubicSpline: fpar_tmp = (Int)fit_param.npiece; break; case BaselineType_kSinusoid: fpar_tmp.resize(fit_param.nwave.size()); fpar_tmp = casacore::Vector<Int>(std::vector<int>(fit_param.nwave.begin(), fit_param.nwave.end())); break; default: throw(AipsError("Unsupported baseline type.")); } if(bltype == BaselineType_kSinusoid){ //Resize the tmp func_param to accomodate number of waves if the size is smaller if (fpar_mtx.shape()[1] < fpar_tmp.size()) { fpar_mtx.resize(IPosition(2, num_pol, fpar_tmp.size()), true); } for(size_t i = 0; i < fpar_tmp.size(); ++i){ fpar_mtx[i][ipol] = fpar_tmp[i]; } } else{ fpar_mtx[0][ipol] = fpar_tmp; } if (num_coeff > num_coeff_max) { num_coeff_max = num_coeff; } SakuraAlignedArray<double> coeff(num_coeff); // CAUTION!!! // data() method must be used with special care!!! double *coeff_data = coeff.data(); string get_coeff_funcname; switch (bltype) { case BaselineType_kPolynomial: case BaselineType_kChebyshev: status = LIBSAKURA_SYMBOL(LSQFitPolynomialFloat)( context, fit_param.order, num_chan, spec_data, mask_data, fit_param.clip_threshold_sigma, fit_param.num_fitting_max, num_coeff, coeff_data, nullptr, nullptr, mask_after_clipping_data, &rms, &bl_status); for (size_t i = 0; i < num_chan; ++i) { if (mask_data[i] == false) { final_mask[ipol] += 1; } if (mask_after_clipping_data[i] == false) { final_mask_after_clipping[ipol] += 1; } } get_coeff_funcname = "sakura_LSQFitPolynomialFloat"; break; case BaselineType_kCubicSpline: status = LIBSAKURA_SYMBOL(LSQFitCubicSplineFloat)( context, fit_param.npiece, num_chan, spec_data, mask_data, fit_param.clip_threshold_sigma, fit_param.num_fitting_max, reinterpret_cast<double (*)[4]>(coeff_data), nullptr, nullptr, mask_after_clipping_data, &rms, boundary_data, &bl_status); for (size_t i = 0; i < num_chan; ++i) { if (mask_data[i] == false) { final_mask[ipol] += 1; } if (mask_after_clipping_data[i] == false) { final_mask_after_clipping[ipol] += 1; } } get_coeff_funcname = "sakura_LSQFitCubicSplineFloat"; break; case BaselineType_kSinusoid: status = LIBSAKURA_SYMBOL(LSQFitSinusoidFloat)( context, fit_param.nwave.size(), &fit_param.nwave[0], num_chan, spec_data, mask_data, fit_param.clip_threshold_sigma, fit_param.num_fitting_max, num_coeff, coeff_data, nullptr, nullptr, mask_after_clipping_data, &rms, &bl_status); for (size_t i = 0; i < num_chan; ++i) { if (mask_data[i] == false) { final_mask[ipol] += 1; } if (mask_after_clipping_data[i] == false) { final_mask_after_clipping[ipol] += 1; } } get_coeff_funcname = "sakura_LSQFitSinusoidFloat"; break; default: throw(AipsError("Unsupported baseline type.")); } check_sakura_status(get_coeff_funcname, status); //For separating of cubic spline in pieces const auto num_ffpar = get_num_coeff_bloutput(fit_param.baseline_type, fit_param.npiece, num_ffpar_max); ffpar_mtx_tmp[ipol].clear(); for (size_t ipiece = 0; ipiece < num_ffpar; ++ipiece) { ffpar_mtx_tmp[ipol].push_back(boundary_data[ipiece]); } coeff_mtx_tmp[ipol].clear(); for (size_t icoeff = 0; icoeff < num_coeff; ++icoeff) { coeff_mtx_tmp[ipol].push_back(coeff_data[icoeff]); } Vector<uInt> masklist; get_masklist_from_mask(num_chan, mask_after_clipping_data, masklist); if (masklist.size() > num_masklist_max) { num_masklist_max = masklist.size(); } masklist_mtx_tmp[ipol].clear(); for (size_t imask = 0; imask < masklist.size(); ++imask) { masklist_mtx_tmp[ipol].push_back(masklist[imask]); } string subtract_funcname; switch (fit_param.baseline_type) { case BaselineType_kPolynomial: case BaselineType_kChebyshev: status = LIBSAKURA_SYMBOL(SubtractPolynomialFloat)( context, num_chan, spec_data, num_coeff, coeff_data, spec_data); subtract_funcname = "sakura_SubtractPolynomialFloat"; break; case BaselineType_kCubicSpline: status = LIBSAKURA_SYMBOL(SubtractCubicSplineFloat)( context, num_chan, spec_data, fit_param.npiece, reinterpret_cast<double (*)[4]>(coeff_data), boundary_data, spec_data); subtract_funcname = "sakura_SubtractCubicSplineFloat"; break; case BaselineType_kSinusoid: status = LIBSAKURA_SYMBOL(SubtractSinusoidFloat)( context, num_chan, spec_data, fit_param.nwave.size(), &fit_param.nwave[0], num_coeff, coeff_data, spec_data); subtract_funcname = "sakura_SubtractSinusoidFloat"; break; default: throw(AipsError("Unsupported baseline type.")); } check_sakura_status(subtract_funcname, status); rms_mtx[0][ipol] = rms; cthres_mtx[0][ipol] = fit_param.clip_threshold_sigma; citer_mtx[0][ipol] = (uInt)fit_param.num_fitting_max - 1; uself_mtx[0][ipol] = (Bool)fit_param.line_finder.use_line_finder; lfthres_mtx[0][ipol] = fit_param.line_finder.threshold; lfavg_mtx[0][ipol] = fit_param.line_finder.chan_avg_limit; for (size_t iedge = 0; iedge < 2; ++iedge) { lfedge_mtx[iedge][ipol] = fit_param.line_finder.edge[iedge]; } } else { string subtract_funcname; switch (fit_param.baseline_type) { case BaselineType_kPolynomial: case BaselineType_kChebyshev: status = LIBSAKURA_SYMBOL(LSQFitPolynomialFloat)( context, fit_param.order, num_chan, spec_data, mask_data, fit_param.clip_threshold_sigma, fit_param.num_fitting_max, num_coeff, nullptr, nullptr, spec_data, mask_after_clipping_data, &rms, &bl_status); subtract_funcname = "sakura_LSQFitPolynomialFloat"; break; case BaselineType_kCubicSpline: status = LIBSAKURA_SYMBOL(LSQFitCubicSplineFloat)( context, fit_param.npiece, num_chan, spec_data, mask_data, fit_param.clip_threshold_sigma, fit_param.num_fitting_max, nullptr, nullptr, spec_data, mask_after_clipping_data, &rms, boundary_data, &bl_status); subtract_funcname = "sakura_LSQFitCubicSplineFloat"; break; case BaselineType_kSinusoid: status = LIBSAKURA_SYMBOL(LSQFitSinusoidFloat)( context, fit_param.nwave.size(), &fit_param.nwave[0], num_chan, spec_data, mask_data, fit_param.clip_threshold_sigma, fit_param.num_fitting_max, num_coeff, nullptr, nullptr, spec_data, mask_after_clipping_data, &rms, &bl_status); subtract_funcname = "sakura_LSQFitSinusoidFloat"; break; default: throw(AipsError("Unsupported baseline type.")); } check_sakura_status(subtract_funcname, status); } // set back a spectrum to data cube if (do_subtract) { set_spectrum_to_cube(data_chunk, irow, ipol, num_chan, spec_data); } if (update_weight) { compute_weight(num_chan, spec_data, mask_data, weight); set_weight_to_matrix(weight_matrix, irow, ipol, weight.at(var_index)); } } // end of polarization loop // output results of fitting if (num_apply_true == 0) continue; Array<Float> ffpar_mtx(IPosition(2, num_pol, num_ffpar_max)); set_matrix_for_bltable<double, Float>(num_pol, num_ffpar_max, ffpar_mtx_tmp, ffpar_mtx); Array<uInt> masklist_mtx(IPosition(2, num_pol, num_masklist_max)); set_matrix_for_bltable<uInt, uInt>(num_pol, num_masklist_max, masklist_mtx_tmp, masklist_mtx); Array<Float> coeff_mtx(IPosition(2, num_pol, num_coeff_max)); set_matrix_for_bltable<double, Float>(num_pol, num_coeff_max, coeff_mtx_tmp, coeff_mtx); Matrix<uInt> masklist_mtx2 = masklist_mtx; Matrix<Bool> apply_mtx2 = apply_mtx; if (write_baseline_table) { bt->appenddata((uInt)scans[irow], (uInt)beams[irow], (uInt)antennas[irow], (uInt)data_spw[irow], 0, times[irow], apply_mtx, bltype_mtx, fpar_mtx, ffpar_mtx, masklist_mtx, coeff_mtx, rms_mtx, (uInt)num_chan, cthres_mtx, citer_mtx, uself_mtx, lfthres_mtx, lfavg_mtx, lfedge_mtx); } if (write_baseline_text) { for (size_t ipol = 0; ipol < num_pol; ++ipol) { if (apply_mtx2(ipol, 0) == false) continue; ofs_txt << "Scan" << '[' << (uInt)scans[irow] << ']' << ' ' << "Beam" << '[' << (uInt)beams[irow] << ']' << ' ' << "Spw" << '[' << (uInt)data_spw[irow] << ']' << ' ' << "Pol" << '[' << ipol << ']' << ' ' << "Time" << '[' << MVTime(times[irow]/ 24. / 3600.).string(MVTime::YMD, 8) << ']' << endl; ofs_txt << endl; ofs_txt << "Fitter range = " << '['; for (size_t imasklist = 0; imasklist < num_masklist_max / 2; ++imasklist) { if (imasklist == 0) { ofs_txt << '[' << masklist_mtx2(ipol, 2 * imasklist) << ';' << masklist_mtx2(ipol, 2 * imasklist + 1) << ']'; } if (imasklist >= 1 && (0 != masklist_mtx2(ipol, 2 * imasklist) && 0 != masklist_mtx2(ipol, 2 * imasklist + 1))) { ofs_txt << ",[" << masklist_mtx2(ipol, 2 * imasklist) << ',' << masklist_mtx2(ipol, 2 * imasklist + 1) << ']'; } } ofs_txt << ']' << endl; ofs_txt << endl; Matrix<uInt> bltype_mtx2 = bltype_mtx[0][ipol]; Matrix<Int> fpar_mtx2 = fpar_mtx[0][ipol]; Matrix<Float> rms_mtx2 = rms_mtx[0][ipol]; string bltype_name; string blparam_name = "order"; uInt bltype_idx = bltype_mtx2(0, 0); if (bltype_idx == BaselineType_kPolynomial) { bltype_name = "poly"; } else if (bltype_idx == BaselineType_kChebyshev) { bltype_name = "chebyshev"; } else if (bltype_idx == BaselineType_kCubicSpline) { bltype_name = "cspline"; blparam_name = "npiece"; } else if (bltype_idx == BaselineType_kSinusoid) { bltype_name = "sinusoid"; blparam_name = "nwave"; } if(bltype_idx == BaselineType_kSinusoid){ ofs_txt << "Baseline parameters Function = " << bltype_name << " " << blparam_name << " = "; Matrix<Int> fpar_mtx3 = fpar_mtx; size_t nwave_size = fpar_mtx3.ncolumn(); for (size_t num = 0; num < nwave_size; ++num) { ofs_txt << fpar_mtx3(ipol, num) << " "; } ofs_txt << endl; } else{ ofs_txt << "Baseline parameters Function = " << bltype_name << " " << blparam_name << " = " << fpar_mtx2(0, 0) << endl; } ofs_txt << endl; ofs_txt << "Results of baseline fit" << endl; ofs_txt << endl; Matrix<Float> coeff_mtx2 = coeff_mtx; for (size_t icoeff = 0; icoeff < numcoeff[ipol]; ++icoeff) { ofs_txt << "p" << icoeff << " = " << setprecision(8) << coeff_mtx2(ipol, icoeff) << " "; } ofs_txt << endl; ofs_txt << endl; ofs_txt << "rms = "; ofs_txt << setprecision(8) << rms_mtx2(0, 0) << endl; ofs_txt << endl; ofs_txt << "Number of clipped channels = " << final_mask_after_clipping[ipol] - final_mask[ipol] << endl; ofs_txt << endl; ofs_txt << "------------------------------------------------------" << endl; ofs_txt << endl; } } if (write_baseline_csv) { for (size_t ipol = 0; ipol < num_pol; ++ipol) { if (apply_mtx2(ipol, 0) == false) continue; ofs_csv << (uInt)scans[irow] << ',' << (uInt)beams[irow] << ',' << (uInt)data_spw[irow] << ',' << ipol << ',' << setprecision(12) << times[irow] << ','; ofs_csv << '['; for (size_t imasklist = 0; imasklist < num_masklist_max / 2; ++imasklist) { if (imasklist == 0) { ofs_csv << '[' << masklist_mtx2(ipol, 2 * imasklist) << ';' << masklist_mtx2(ipol, 2 * imasklist + 1) << ']'; } if (imasklist >= 1 && (0 != masklist_mtx2(ipol, 2 * imasklist) && 0 != masklist_mtx2(ipol, 2 * imasklist + 1))) { ofs_csv << ";[" << masklist_mtx2(ipol, 2 * imasklist) << ';' << masklist_mtx2(ipol, 2 * imasklist + 1) << ']'; } } ofs_csv << ']' << ','; Matrix<uInt> bltype_mtx2 = bltype_mtx[0][ipol]; string bltype_name; uInt bltype_idx = bltype_mtx2(0, 0); if (bltype_idx == BaselineType_kPolynomial) { bltype_name = "poly"; } else if (bltype_idx == BaselineType_kChebyshev) { bltype_name = "chebyshev"; } else if (bltype_idx == BaselineType_kCubicSpline) { bltype_name = "cspline"; } else if (bltype_idx== BaselineType_kSinusoid) { bltype_name = "sinusoid"; } Matrix<Int> fpar_mtx2 = fpar_mtx; Matrix<Float> coeff_mtx2 = coeff_mtx; if (bltype_idx == BaselineType_kSinusoid) { size_t nwave_size = fpar_mtx2.ncolumn(); ofs_csv << bltype_name.c_str() << ',' << fpar_mtx2(ipol, 0); for (size_t i = 1; i < nwave_size; ++i) { ofs_csv << ';' << fpar_mtx2(ipol, i); } ofs_csv << ','; } else { ofs_csv << bltype_name.c_str() << ',' << fpar_mtx2(ipol, 0) << ','; } for (size_t icoeff = 0; icoeff < numcoeff[ipol]; ++icoeff) { ofs_csv << setprecision(8) << coeff_mtx2(ipol, icoeff) << ','; } Matrix<Float> rms_mtx2 = rms_mtx; ofs_csv << setprecision(8) << rms_mtx2(ipol, 0) << ','; ofs_csv << final_mask_after_clipping[ipol] - final_mask[ipol]; ofs_csv << endl; } } } // end of chunk row loop // write back data and flag cube to VisBuffer if (update_weight) { sdh_->fillCubeToOutputMs(vb, data_chunk, &flag_chunk, &weight_matrix); } else { sdh_->fillCubeToOutputMs(vb, data_chunk, &flag_chunk); } } // end of vi loop } // end of chunk loop if (write_baseline_csv) { ofs_csv.close(); } if (write_baseline_text) { ofs_txt.close(); } if (write_baseline_table) { bt->save(bloutputname_table); delete bt; } finalize_process(); // destroy baseline contexts map<size_t const, std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> >::iterator ctxiter = context_reservoir.begin(); while (ctxiter != context_reservoir.end()) { destroy_baseline_contexts (context_reservoir[(*ctxiter).first]); ++ctxiter; } } //end subtractBaselineVariable list<pair<size_t, size_t>> SingleDishMS::findLineAndGetRanges(size_t const num_data, float const* data, bool * mask, float const threshold, int const avg_limit, int const minwidth, vector<int> const& edge, bool const invert) { // input value check AlwaysAssert(minwidth > 0, AipsError); AlwaysAssert(avg_limit >= 0, AipsError); size_t max_iteration = 10; size_t maxwidth = num_data; AlwaysAssert(maxwidth > static_cast<size_t>(minwidth), AipsError); // edge handling pair<size_t, size_t> lf_edge; if (edge.size() == 0) { lf_edge = pair<size_t, size_t>(0, 0); } else if (edge.size() == 1) { AlwaysAssert(edge[0] >= 0, AipsError); lf_edge = pair<size_t, size_t>(static_cast<size_t>(edge[0]), static_cast<size_t>(edge[0])); } else { AlwaysAssert(edge[0] >= 0 && edge[1] >= 0, AipsError); lf_edge = pair<size_t, size_t>(static_cast<size_t>(edge[0]), static_cast<size_t>(edge[1])); } // line detection list<pair<size_t, size_t>> line_ranges = linefinder::MADLineFinder(num_data, data, mask, threshold, max_iteration, static_cast<size_t>(minwidth), maxwidth, static_cast<size_t>(avg_limit), lf_edge); // debug output LogIO os(_ORIGIN); os << LogIO::DEBUG1 << line_ranges.size() << " lines found: "; for (list<pair<size_t, size_t>>::iterator iter = line_ranges.begin(); iter != line_ranges.end(); ++iter) { os << "[" << (*iter).first << ", " << (*iter).second << "] "; } os << LogIO::POST; if (invert) { // eliminate edge channels from output mask if (lf_edge.first > 0) line_ranges.push_front(pair<size_t, size_t>(0, lf_edge.first - 1)); if (lf_edge.second > 0) line_ranges.push_back( pair<size_t, size_t>(num_data - lf_edge.second, num_data - 1)); } return line_ranges; } void SingleDishMS::findLineAndGetMask(size_t const num_data, float const* data, bool const* in_mask, float const threshold, int const avg_limit, int const minwidth, vector<int> const& edge, bool const invert, bool* out_mask) { // copy input mask to output mask vector if necessary if (in_mask != out_mask) { for (size_t i = 0; i < num_data; ++i) { out_mask[i] = in_mask[i]; } } // line finding list<pair<size_t, size_t>> line_ranges = findLineAndGetRanges(num_data, data, out_mask, threshold, avg_limit, minwidth, edge, invert); // line mask creation (do not initialize in case of baseline mask) linefinder::getMask(num_data, out_mask, line_ranges, invert, !invert); } void SingleDishMS::smooth(string const &kernelType, float const kernelWidth, string const &columnName, string const &outMSName) { LogIO os(_ORIGIN); os << "Input parameter summary:" << endl << " kernelType = " << kernelType << endl << " kernelWidth = " << kernelWidth << endl << " columnName = " << columnName << endl << " outMSName = " << outMSName << LogIO::POST; // Initialization doSmoothing_ = true; prepare_for_process(columnName, outMSName); // configure smoothing sdh_->setSmoothing(kernelType, kernelWidth); sdh_->initializeSmoothing(); // get VI/VB2 access vi::VisibilityIterator2 *visIter = sdh_->getVisIter(); visIter->setRowBlocking(kNRowBlocking); vi::VisBuffer2 *vb = visIter->getVisBuffer(); double startTime = gettimeofday_sec(); for (visIter->originChunks(); visIter->moreChunks(); visIter->nextChunk()) { for (visIter->origin(); visIter->more(); visIter->next()) { sdh_->fillOutputMs(vb); } } double endTime = gettimeofday_sec(); os << LogIO::DEBUGGING << "Elapsed time for VI/VB loop: " << endTime - startTime << " sec" << LogIO::POST; // Finalization finalize_process(); } void SingleDishMS::atmcor(Record const &config, string const &columnName, string const &outMSName) { LogIO os(_ORIGIN); os << LogIO::DEBUGGING << "Input parameter summary:" << endl << " columnName = " << columnName << endl << " outMSName = " << outMSName << LogIO::POST; // Initialization doAtmCor_ = true; atmCorConfig_ = config; os << LogIO::DEBUGGING << "config summry:"; atmCorConfig_.print(os.output(), 25, " "); os << LogIO::POST; Block<Int> sortCols(4); sortCols[0] = MS::OBSERVATION_ID; sortCols[1] = MS::ARRAY_ID; sortCols[2] = MS::FEED1; sortCols[3] = MS::DATA_DESC_ID; prepare_for_process(columnName, outMSName, sortCols, False); // get VI/VB2 access vi::VisibilityIterator2 *visIter = sdh_->getVisIter(); // for parallel processing: set row blocking (common multiple of 3 and 4) // TODO: optimize row blocking size constexpr rownr_t kNrowBlocking = 360u; std::vector<Int> antenna1 = ScalarColumn<Int>(visIter->ms(), "ANTENNA1").getColumn().tovector(); std::sort(antenna1.begin(), antenna1.end()); auto const result = std::unique(antenna1.begin(), antenna1.end()); Int const nAntennas = std::distance(antenna1.begin(), result); visIter->setRowBlocking(kNrowBlocking * nAntennas); os << "There are " << nAntennas << " antennas in MAIN table. " << "Set row-blocking size " << kNrowBlocking * nAntennas << LogIO::POST; vi::VisBuffer2 *vb = visIter->getVisBuffer(); double startTime = gettimeofday_sec(); for (visIter->originChunks(); visIter->moreChunks(); visIter->nextChunk()) { for (visIter->origin(); visIter->more(); visIter->next()) { sdh_->fillOutputMs(vb); } } double endTime = gettimeofday_sec(); os << LogIO::DEBUGGING << "Elapsed time for VI/VB loop: " << endTime - startTime << " sec" << LogIO::POST; // Finalization finalize_process(); } bool SingleDishMS::importAsap(string const &infile, string const &outfile, bool const parallel) { bool status = true; try { SingleDishMSFiller<Scantable2MSReader> filler(infile, parallel); filler.fill(); filler.save(outfile); } catch (AipsError &e) { LogIO os(_ORIGIN); os << LogIO::SEVERE << "Exception occurred." << LogIO::POST; os << LogIO::SEVERE << "Original Message: \n" << e.getMesg() << LogIO::POST; os << LogIO::DEBUGGING << "Detailed Stack Trace: \n" << e.getStackTrace() << LogIO::POST; status = false; } catch (...) { LogIO os(_ORIGIN); os << LogIO::SEVERE << "Unknown exception occurred." << LogIO::POST; status = false; } return status; } bool SingleDishMS::importNRO(string const &infile, string const &outfile, bool const parallel) { bool status = true; try { SingleDishMSFiller<NRO2MSReader> filler(infile, parallel); filler.fill(); filler.save(outfile); } catch (AipsError &e) { LogIO os(_ORIGIN); os << LogIO::SEVERE << "Exception occurred." << LogIO::POST; os << LogIO::SEVERE << "Original Message: \n" << e.getMesg() << LogIO::POST; os << LogIO::DEBUGGING << "Detailed Stack Trace: \n" << e.getStackTrace() << LogIO::POST; status = false; } catch (...) { LogIO os(_ORIGIN); os << LogIO::SEVERE << "Unknown exception occurred." << LogIO::POST; status = false; } return status; } } // End of casa namespace.